Repository: cogent3/cogent3
Pull request #2141 (2024-11-29 01:49:11)

ENH: major performance improvement in pairwise distance calculations
## Summary by Sourcery

Refactor the test suite to use pytest fixtures and parameterized tests, enhancing modularity and readability. Implement numba-accelerated functions for pairwise distance calculations, significantly improving performance for various distance metrics including JC69, TN93, paralinear, hamming, and proportion difference.

Enhancements:
- Refactor test suite to use pytest fixtures for better modularity and reusability.
- Optimize pairwise distance calculations by implementing numba-accelerated functions for various distance metrics including JC69, TN93, paralinear, hamming, and proportion difference.

Tests:
- Convert existing test cases to use pytest framework, improving test structure and readability.
- Add parameterized tests for numba-accelerated distance calculations to ensure accuracy and performance.


Patch:
diff --git a/src/cogent3/core/new_alignment.py b/src/cogent3/core/new_alignment.py
index 7cf562157..4aac71185 100644
--- a/src/cogent3/core/new_alignment.py
+++ b/src/cogent3/core/new_alignment.py
@@ -5008,39 +5008,35 @@ class Alignment(SequenceCollection):
     def distance_matrix(
         self,
         calc: str = "pdist",
-        show_progress: bool = False,
         drop_invalid: bool = False,
+        parallel: bool = False,
     ):
         """Returns pairwise distances between sequences.
 
         Parameters
         ----------
         calc
-            a pairwise distance calculator or name of one. For options see
-            cogent3.evolve.fast_distance.available_distances
-        show_progress
-            controls progress display for distance calculation
+            a pairwise distance calculator name. Presently only
+            'pdist', 'jc69', 'tn93', 'hamming', 'paralinear' are supported.
         drop_invalid
             If True, sequences for which a pairwise distance could not be
             calculated are excluded. If False, an ArithmeticError is raised if
             a distance could not be computed on observed data.
         """
-        from cogent3.evolve.fast_distance import get_distance_calculator
+        from cogent3.evolve.pairwise_distance_numba import get_distance_calculator
 
+        calc = get_distance_calculator(
+            calc,
+        )
         try:
-            calculator = get_distance_calculator(
-                calc,
-                moltype=self.moltype,
-                alignment=self,
-                invalid_raises=not drop_invalid,
-            )
-            calculator.run(show_progress=show_progress)
+            result = calc(self, invalid_raises=not drop_invalid, parallel=parallel)
         except ArithmeticError as e:
             msg = "not all pairwise distances could be computed, try drop_invalid=True"
             raise ArithmeticError(msg) from e
-        result = calculator.get_pairwise_distances()
+
         if drop_invalid:
             result = result.drop_invalid()
+
         return result
 
     @UI.display_wrap
@@ -5049,6 +5045,7 @@ class Alignment(SequenceCollection):
         calc: str = "pdist",
         bootstrap: OptInt = None,
         drop_invalid: bool = False,
+        parallel: bool = False,
         show_progress: bool = False,
         ui=None,  # refactor: type hint
     ):
@@ -5067,6 +5064,8 @@ class Alignment(SequenceCollection):
             If True, sequences for which a pairwise distance could not be
             calculated are excluded. If False, an ArithmeticError is raised if
             a distance could not be computed on observed data.
+        parallel
+            parallel execution of distance calculations
         show_progress
             controls progress display for distance calculation
 
@@ -5086,7 +5085,7 @@ class Alignment(SequenceCollection):
         from cogent3.phylo.nj import gnj
 
         dm = self.distance_matrix(
-            calc=calc, show_progress=show_progress, drop_invalid=drop_invalid
+            calc=calc, drop_invalid=drop_invalid, parallel=parallel
         )
         results = gnj(dm, keep=1, show_progress=show_progress)
         kept_names = set(dm.template.names[0])
@@ -5098,7 +5097,7 @@ class Alignment(SequenceCollection):
                 b = subaln.sample(with_replacement=True)
                 try:
                     bdist = b.distance_matrix(
-                        calc=calc, show_progress=show_progress, drop_invalid=False
+                        calc=calc, drop_invalid=drop_invalid, parallel=parallel
                     )
                 except ArithmeticError:
                     # all sequences must have a pairwise distance to allow
diff --git a/src/cogent3/core/new_moltype.py b/src/cogent3/core/new_moltype.py
index bff1cf40b..e1c70d871 100644
--- a/src/cogent3/core/new_moltype.py
+++ b/src/cogent3/core/new_moltype.py
@@ -617,6 +617,8 @@ class MolType:
             raise new_alphabet.AlphabetError(
                 f"{seq[:4]!r} not valid for moltype {self.name!r}"
             )
+        if not self.is_nucleic:
+            raise MolTypeError(f"{self.name!r} cannot complement")
         return self._complement(seq)
 
     @complement.register
diff --git a/src/cogent3/evolve/fast_distance.py b/src/cogent3/evolve/fast_distance.py
index b49bb1e9e..ce57b9b5b 100644
--- a/src/cogent3/evolve/fast_distance.py
+++ b/src/cogent3/evolve/fast_distance.py
@@ -2,6 +2,7 @@ from collections import defaultdict, namedtuple
 from numbers import Number
 from typing import Tuple
 
+import numba
 import numpy
 from numpy import array, diag, dot, eye, float64, int32, log, sqrt, zeros
 from numpy.linalg import det, inv
@@ -13,9 +14,19 @@ from cogent3.util.dict_array import DictArray
 from cogent3.util.misc import get_object_provenance
 from cogent3.util.progress_display import display_wrap
 
-from .pairwise_distance_numba import fill_diversity_matrix
 
-# pending addition of protein distance metrics
+@numba.jit(cache=True)
+def fill_diversity_matrix(matrix, seq1, seq2):  # pragma: no cover
+    """fills the diversity matrix for valid positions.
+
+    Assumes the provided sequences have been converted to indices with
+    invalid characters being negative numbers (use get_moltype_index_array
+    plus seq_to_indices)."""
+
+    for i in range(len(seq1)):
+        if seq1[i] < 0 or seq2[i] < 0:
+            continue
+        matrix[seq1[i], seq2[i]] += 1.0
 
 
 def _same_moltype(ref, query):
@@ -25,22 +36,19 @@ def _same_moltype(ref, query):
 
 def get_pyrimidine_indices(moltype):
     """returns pyrimidine indices for the moltype"""
-    states = list(moltype)
-    if _same_moltype(RNA, moltype):
-        return list(map(states.index, "CU"))
-    elif _same_moltype(DNA, moltype):
-        return list(map(states.index, "CT"))
-    else:
+    alpha = moltype.alphabet
+    if len(alpha) != 4:
         raise RuntimeError("Non-nucleic acid MolType")
+    pyrimidines = "CT" if moltype.label == "dna" else "CU"
+    return alpha.to_indices(pyrimidines).tolist()
 
 
 def get_purine_indices(moltype):
     """returns purine indices for the moltype"""
-    states = list(moltype)
-    if not _same_moltype(RNA, moltype) and not _same_moltype(DNA, moltype):
+    alpha = moltype.alphabet
+    if len(alpha) != 4:
         raise RuntimeError("Non-nucleic acid MolType")
-
-    return list(map(states.index, "AG"))
+    return alpha.to_indices("AG").tolist()
 
 
 def get_matrix_diff_coords(indices):
@@ -80,18 +88,6 @@ def seq_to_indices(seq, char_to_index):
     return char_to_index.take(ords)
 
 
-def _fill_diversity_matrix(matrix, seq1, seq2):
-    """fills the diversity matrix for valid positions.
-
-    Assumes the provided sequences have been converted to indices with
-    invalid characters being negative numbers (use get_moltype_index_array
-    plus seq_to_indices)."""
-    paired = array([seq1, seq2]).T
-    paired = paired[paired.min(axis=1) >= 0]
-    for i in range(len(paired)):
-        matrix[paired[i][0], paired[i][1]] += 1
-
-
 def _hamming(matrix):
     """computes the edit distance
     Parameters
diff --git a/src/cogent3/evolve/pairwise_distance_numba.py b/src/cogent3/evolve/pairwise_distance_numba.py
index ae88d74e9..e417aaf2b 100644
--- a/src/cogent3/evolve/pairwise_distance_numba.py
+++ b/src/cogent3/evolve/pairwise_distance_numba.py
@@ -1,18 +1,526 @@
-from numba import njit
+import itertools
+import typing
 
-# turn off code coverage as njit-ted code not accessible to coverage
+import numba
+import numba.types as numba_types
+import numpy
+
+from cogent3.core import new_moltype
+from cogent3.evolve.fast_distance import DistanceMatrix
+
+if typing.TYPE_CHECKING:
+    # Forward reference to avoid circular import
+    from cogent3.core.new_alignment import Alignment
+
+# turn off code coverage as jit-ted code not accessible to coverage
 
 
 # fills in a diversity matrix from sequences of integers
-@njit(cache=True)
-def fill_diversity_matrix(matrix, seq1, seq2):  # pragma: no cover
-    """fills the diversity matrix for valid positions.
 
-    Assumes the provided sequences have been converted to indices with
-    invalid characters being negative numbers (use get_moltype_index_array
-    plus seq_to_indices)."""
 
+def _is_nucleic(moltype):
+    # todo delete when new_type alignments and moltypes are the norm
+    try:
+        _ = moltype.complement("A")
+    except (ValueError, new_moltype.MolTypeError):
+        return False
+    return True
+
+
+@numba.jit
+def num_diffs_and_valid(
+    seq1: numba_types.uint8[:],
+    seq2: numba_types.uint8[:],
+    num_states: numba_types.uint8,
+) -> numba_types.int32[:]:  # pragma: no cover
+    """returns the number of differences and the number of valid positions"""
+    num_valid = 0
+    num_diff = 0
     for i in range(len(seq1)):
-        if seq1[i] < 0 or seq2[i] < 0:
+        if seq1[i] >= num_states or seq2[i] >= num_states:
             continue
-        matrix[seq1[i], seq2[i]] += 1.0
+        num_valid += 1
+        if seq1[i] != seq2[i]:
+            num_diff += 1
+
+    return numpy.array([num_diff, num_valid], dtype=numpy.int32)
+
+
+@numba.jit(cache=True)
+def _get_matrix_and_counts(matrix, seq1, seq2):  # pragma: no cover
+    """fills the diversity matrix for valid positions.
+
+    Parameters
+    ----------
+    matrix
+        Counts of states
+    seq1, seq2
+        Sequences to count variations between
+    """
+    num_valid = 0
+    num_diff = 0
+    num_states = matrix.shape[0]
+    for i in range(len(seq1)):
+        if seq1[i] >= num_states or seq2[i] >= num_states:
+            continue
+        num_valid += 1
+        if seq1[i] != seq2[i]:
+            num_diff += 1
+        matrix[seq1[i], seq2[i]] += 1
+    return matrix, num_diff, num_valid
+
+
+@numba.jit(parallel=True)
+def jc69_dist_matrix(array_seqs, num_states, parallel=True):  # pragma: no cover
+    """Compute JC69 distances between all sequence pairs.
+
+    Parameters
+    ----------
+    array_seqs
+        2D array of sequences
+    num_states
+        Number of states, e.g. 4 for nucleotides
+    parallel
+        Runs in parallel if True, otherwise reduces to single thread
+
+    Returns
+    -------
+    array
+        Matrix of pairwise JC69 distances
+
+    Notes
+    -----
+    If the proportion of differences is >= 0.75, the distance is set
+    to numpy.nan.
+    """
+    num_threads = numba.get_num_threads()
+    if not parallel:
+        # numba still gets all the threads it can, but it only
+        # uses the number specified by set_num_threads
+        # we restore the original number of threads at the end
+        numba.set_num_threads(1)
+
+    n_seqs = array_seqs.shape[0]
+    dists = numpy.zeros((n_seqs, n_seqs), dtype=numpy.float32)
+    nan = numpy.float32(numpy.nan)
+    zero = numpy.float32(0.0)
+    frac = numpy.float32(num_states / (num_states - 1))
+    dist_scale = numpy.float32((num_states - 1) / -num_states)
+
+    # outer loop is parallelised
+    for i in numba.prange(n_seqs - 1):
+        for j in range(i + 1, n_seqs):
+            num_diffs, num_valid = num_diffs_and_valid(
+                array_seqs[i], array_seqs[j], num_states
+            )
+
+            if num_valid == 0:
+                d = nan
+            elif num_diffs == 0:
+                d = zero
+            else:
+                p = num_diffs / num_valid
+                d = nan if p >= 0.75 else dist_scale * numpy.log(1.0 - frac * p)
+            dists[i, j] = dists[j, i] = d
+
+    if not parallel:
+        # restore original number of threads
+        numba.set_num_threads(num_threads)
+    return dists
+
+
+def jc69(
+    aln: "Alignment", invalid_raises: bool = False, parallel: bool = False
+) -> DistanceMatrix:
+    """returns JC69 pairwise distances for an alignment
+
+    Parameters
+    ----------
+    aln
+        Alignment instance
+    invalid_raises
+        If True and nan in result, raises an ArithmeticError
+    parallel
+        If True, uses parallel processing via numba.
+    """
+
+    num_states = len(aln.moltype.alphabet)
+    mat = jc69_dist_matrix(aln.array_seqs, num_states, parallel=parallel)
+    if invalid_raises and numpy.isnan(mat).any():
+        raise ArithmeticError("nan's in matrix")
+    return DistanceMatrix.from_array_names(mat, aln.names)
+
+
+@numba.jit
+def _calc_tn93_dist(
+    matrix,
+    pur_coords,
+    pyr_coords,
+    tv_coords,
+    freq_purs,
+    freq_pyrs,
+    coeff1,
+    coeff2,
+    coeff3,
+    total,
+):  # pragma: no cover
+    pur_ts_diffs = matrix.take(pur_coords).sum() / total
+    pyr_ts_diffs = matrix.take(pyr_coords).sum() / total
+    tv_diffs = matrix.flatten().take(tv_coords).sum() / total
+    term1 = 1 - pur_ts_diffs / coeff1 - tv_diffs / (2 * freq_purs)
+    term2 = 1 - pyr_ts_diffs / coeff2 - tv_diffs / (2 * freq_pyrs)
+    term3 = 1 - tv_diffs / (2 * freq_purs * freq_pyrs)
+
+    if term1 <= 0 or term2 <= 0 or term3 <= 0:  # log will fail
+        return numpy.float32(numpy.nan)
+
+    return (
+        -coeff1 * numpy.log(term1)
+        - coeff2 * numpy.log(term2)
+        - coeff3 * numpy.log(term3)
+    )
+
+
+@numba.jit(parallel=True)
+def tn93_dist_matrix(
+    array_seqs,
+    num_states,
+    pur_freqs,
+    pyr_freqs,
+    pur_coords,
+    pyr_coords,
+    tv_coords,
+    coeff1,
+    coeff2,
+    coeff3,
+    parallel=True,
+):  # pragma: no cover
+    num_threads = numba.get_num_threads()
+    if not parallel:
+        numba.set_num_threads(1)
+    n_seqs = array_seqs.shape[0]
+    dists = numpy.zeros((n_seqs, n_seqs), dtype=numpy.float32)
+    nan = numpy.float32(numpy.nan)
+    zero = numpy.float32(0.0)
+
+    # making working matrices for each thread
+    matrices = numpy.zeros(
+        (numba.get_num_threads(), num_states, num_states), dtype=numpy.int32
+    )
+    # outer parallel loop
+    for i in numba.prange(n_seqs - 1):
+        # get a working matrix for this thread
+        div_matrix = matrices[numba.get_thread_id()]
+        for j in range(i + 1, n_seqs):
+            div_matrix.fill(0)
+            div_matrix, num_diffs, num_valid = _get_matrix_and_counts(
+                div_matrix, array_seqs[i], array_seqs[j]
+            )
+            if num_valid == 0:
+                d = nan
+            elif num_diffs == 0:
+                d = zero
+            else:
+                d = _calc_tn93_dist(
+                    div_matrix.flatten(),
+                    pur_coords,
+                    pyr_coords,
+                    tv_coords,
+                    pur_freqs,
+                    pyr_freqs,
+                    coeff1,
+                    coeff2,
+                    coeff3,
+                    num_valid,
+                )
+            dists[i, j] = dists[j, i] = d
+
+    if not parallel:
+        # restore original number of threads
+        numba.set_num_threads(num_threads)
+    return dists
+
+
+@numba.jit(parallel=True)
+def _count_states(array_seqs, num_states, parallel=True):
+    # temp function to remove when the new Alignment.get_motif_probs()
+    # is refactored to use numba
+    num_threads = numba.get_num_threads()
+    if not parallel:
+        numba.set_num_threads(1)
+    n_seqs = array_seqs.shape[0]
+    # making working matrices for each thread
+    state_counts = numpy.zeros((numba.get_num_threads(), num_states), dtype=numpy.int32)
+    for i in numba.prange(n_seqs):
+        states = state_counts[numba.get_thread_id()]
+        for j in range(len(array_seqs[i])):
+            state = array_seqs[i, j]
+            if state < num_states:
+                states[state] += 1
+
+    if not parallel:
+        # restore original number of threads
+        numba.set_num_threads(num_threads)
+
+    return state_counts.sum(axis=0)
+
+
+def _get_symmetric_within(indices: numpy.ndarray) -> numpy.ndarray:
+    dims = 4, 4
+    coords = numpy.array(
+        [(i, j) for i, j in itertools.product(indices, indices) if i != j]
+    )
+    return numpy.ravel_multi_index(coords.T, dims=dims)
+
+
+def _get_symmetric_between(
+    indices1: numpy.ndarray, indices2: numpy.ndarray
+) -> numpy.ndarray:
+    # coordinates
+    dims = 4, 4
+    coords = numpy.array(
+        list(
+            itertools.chain(
+                itertools.product(indices1, indices2),
+                itertools.product(indices2, indices1),
+            )
+        )
+    )
+    return numpy.ravel_multi_index(coords.T, dims=dims)
+
+
+def tn93(
+    aln: "Alignment", invalid_raises: bool = False, parallel: bool = False
+) -> DistanceMatrix:
+    """returns TN93 pairwise distances for an alignment
+
+    Parameters
+    ----------
+    aln
+        Alignment instance
+    invalid_raises
+        If True and nan in result, raises an ArithmeticError
+    parallel
+        If True, uses parallel processing via numba.
+    """
+    if not _is_nucleic(aln.moltype):
+        raise new_moltype.MolTypeError(
+            f"tn93 distance only works with nucleotide alignments, not {aln.moltype}"
+        )
+
+    alpha = aln.moltype.alphabet
+
+    counts = _count_states(aln.array_seqs, len(alpha), parallel=parallel)
+    freqs = (counts / counts.sum()).astype(numpy.float32)
+
+    pur_indices = alpha.to_indices("AG")
+    pyr_indices = alpha.to_indices("CT")
+
+    # these constants are used across all the pairwise distance calculations
+    pur_freqs = freqs.take(pur_indices).sum()
+    pur_prods = freqs.take(pur_indices).prod()
+    pyr_freqs = freqs.take(pyr_indices).sum()
+    pyr_prods = freqs.take(pyr_indices).prod()
+    coeff1 = 2 * pur_prods / pur_freqs
+    coeff2 = 2 * pyr_prods / pyr_freqs
+    coeff3 = 2 * (
+        pur_freqs * pyr_freqs
+        - (pur_prods * pyr_freqs / pur_freqs)
+        - (pyr_prods * pur_freqs / pyr_freqs)
+    )
+
+    # coords are for the flattened matrix
+    pur_coords = _get_symmetric_within(pur_indices)
+    pyr_coords = _get_symmetric_within(pyr_indices)
+    tv_coords = _get_symmetric_between(pur_indices, pyr_indices)
+
+    mat = tn93_dist_matrix(
+        aln.array_seqs,
+        4,
+        pur_freqs,
+        pyr_freqs,
+        pur_coords,
+        pyr_coords,
+        tv_coords,
+        coeff1,
+        coeff2,
+        coeff3,
+        parallel=parallel,
+    )
+    if invalid_raises and numpy.isnan(mat).any():
+        raise ArithmeticError("nan's in matrix")
+    return DistanceMatrix.from_array_names(mat, aln.names)
+
+
+@numba.jit
+def _paralinear(matrix, num_states):
+    # we replace the missing diagonal states with a
+    # pseudocount of 0.5 then normalise
+    frequency = matrix.astype(numpy.float32)
+    for i in range(num_states):
+        if frequency[i, i] == 0:
+            frequency[i, i] = 0.5
+
+    frequency /= frequency.sum()
+    det = numpy.linalg.det(frequency)
+    if det <= 0:
+        return numpy.float32(numpy.nan)
+
+    # the inverse matrix of frequency, every element is squared
+    denom = numpy.sqrt((frequency.sum(axis=0) * frequency.sum(axis=1)).prod())
+
+    return -numpy.log(det / denom) / num_states
+
+
+@numba.jit(parallel=True)
+def paralinear_distance_matrix(
+    array_seqs,
+    num_states,
+    parallel=True,
+):  # pragma: no cover
+    num_threads = numba.get_num_threads()
+    if not parallel:
+        numba.set_num_threads(1)
+    n_seqs = array_seqs.shape[0]
+    dists = numpy.zeros((n_seqs, n_seqs), dtype=numpy.float32)
+    nan = numpy.float32(numpy.nan)
+    zero = numpy.float32(0.0)
+
+    # making working matrices for each thread
+    matrices = numpy.zeros(
+        (numba.get_num_threads(), num_states, num_states), dtype=numpy.int32
+    )
+    # outer parallel loop
+    for i in numba.prange(n_seqs - 1):
+        # get a working matrix for this thread
+        div_matrix = matrices[numba.get_thread_id()]
+        for j in range(i + 1, n_seqs):
+            div_matrix.fill(0)
+            div_matrix, num_diffs, num_valid = _get_matrix_and_counts(
+                div_matrix, array_seqs[i], array_seqs[j]
+            )
+            if num_valid == 0:
+                d = nan
+            elif num_diffs == 0:
+                d = zero
+            else:
+                d = _paralinear(
+                    div_matrix,
+                    num_states,
+                )
+            dists[i, j] = dists[j, i] = d
+
+    if not parallel:
+        # restore original number of threads
+        numba.set_num_threads(num_threads)
+    return dists
+
+
+def paralinear(
+    aln: "Alignment", invalid_raises: bool = False, parallel: bool = False
+) -> DistanceMatrix:
+    """returns matrix of pairwise paralinear distances
+
+    Parameters
+    ----------
+    aln
+        Alignment instance
+    invalid_raises
+        If True and nan in result, raises an ArithmeticError
+    parallel
+        If True, uses parallel processing via numba.
+
+    Notes
+    -----
+    This is limited to 4-state alphabets for now.
+    """
+    if not _is_nucleic(aln.moltype):
+        raise new_moltype.MolTypeError(
+            f"paralinear distance only works with nucleotide alignments, not {aln.moltype}"
+        )
+
+    num_states = len(aln.moltype.alphabet)
+
+    mat = paralinear_distance_matrix(aln.array_seqs, num_states, parallel=parallel)
+    if invalid_raises and numpy.isnan(mat).any():
+        raise ArithmeticError("nan's in matrix")
+    return DistanceMatrix.from_array_names(mat, aln.names)
+
+
+@numba.jit(parallel=True)
+def simple_distance_matrix(
+    array_seqs, num_states, parallel=True, hamming=True
+):  # pragma: no cover
+    num_threads = numba.get_num_threads()
+    if not parallel:
+        numba.set_num_threads(1)
+
+    n_seqs = array_seqs.shape[0]
+    dists = numpy.zeros((n_seqs, n_seqs), dtype=numpy.float32)
+    # making working matrices for each thread
+    # outer parallel loop
+    for i in numba.prange(n_seqs - 1):
+        for j in range(i + 1, n_seqs):
+            d, num_valid = num_diffs_and_valid(array_seqs[i], array_seqs[j], num_states)
+            if not hamming:
+                d /= num_valid
+            dists[i, j] = dists[j, i] = d
+
+    if not parallel:
+        # restore original number of threads
+        numba.set_num_threads(num_threads)
+    return dists
+
+
+def hamming(aln: "Alignment", parallel: bool = False, **kwargs) -> DistanceMatrix:
+    """returns matrix of hamming distances
+
+    Parameters
+    ----------
+    aln
+        Alignment instance
+    parallel
+        If True, uses parallel processing via numba.
+    """
+    num_states = len(aln.moltype.alphabet)
+    mat = simple_distance_matrix(
+        aln.array_seqs, num_states, parallel=parallel, hamming=True
+    )
+    return DistanceMatrix.from_array_names(mat, aln.names)
+
+
+def pdist(aln: "Alignment", parallel: bool = False, **kwargs) -> DistanceMatrix:
+    """returns matrix of pairwise proportion difference distances
+
+    Parameters
+    ----------
+    aln
+        Alignment instance
+    parallel
+        If True, uses parallel processing via numba.
+    """
+    num_states = len(aln.moltype.alphabet)
+    mat = simple_distance_matrix(
+        aln.array_seqs, num_states, parallel=parallel, hamming=False
+    )
+    return DistanceMatrix.from_array_names(mat, aln.names)
+
+
+_calculators = {
+    "paralinear": paralinear,
+    "jc69": jc69,
+    "tn93": tn93,
+    "hamming": hamming,
+    "pdist": pdist,
+}
+
+
+def get_distance_calculator(name):
+    """returns a pairwise distance calculator
+
+    name is converted to lower case"""
+    name = name.lower()
+    if name not in _calculators:
+        raise ValueError(f'Unknown pairwise distance calculator "{name}"')
+
+    return _calculators[name]


Test patch:
diff --git a/tests/test_core/test_new_alignment.py b/tests/test_core/test_new_alignment.py
index a5e56ba34..015e15f85 100644
--- a/tests/test_core/test_new_alignment.py
+++ b/tests/test_core/test_new_alignment.py
@@ -4246,7 +4246,7 @@ def test_alignment_distance_matrix():
     """Alignment distance_matrix should produce correct scores"""
     data = dict([("s1", "ACGTACGTA"), ("s2", "GTGTACGTA")])
     aln = new_alignment.make_aligned_seqs(data, moltype="dna")
-    dists = aln.distance_matrix(calc="hamming", show_progress=False)
+    dists = aln.distance_matrix(calc="hamming")
     assert dists == {("s1", "s2"): 2.0, ("s2", "s1"): 2.0}
     # and for protein
     aa = aln.get_translation()
diff --git a/tests/test_evolve/test_distance.py b/tests/test_evolve/test_distance.py
index b4e9a51c0..53e41ae33 100644
--- a/tests/test_evolve/test_distance.py
+++ b/tests/test_evolve/test_distance.py
@@ -1,6 +1,4 @@
-import os
 import warnings
-from unittest import TestCase
 
 import numpy
 import pytest
@@ -14,6 +12,7 @@ from cogent3 import (
     make_aligned_seqs,
     make_unaligned_seqs,
 )
+from cogent3.evolve import pairwise_distance_numba as pdist_numba
 from cogent3.evolve.distance import EstimateDistances
 from cogent3.evolve.fast_distance import (
     DistanceMatrix,
@@ -24,18 +23,15 @@ from cogent3.evolve.fast_distance import (
     ProportionIdenticalPair,
     TN93Pair,
     _calculators,
-    _fill_diversity_matrix,
     _hamming,
     _jc69_from_matrix,
     available_distances,
+    fill_diversity_matrix,
     get_distance_calculator,
     get_moltype_index_array,
     seq_to_indices,
 )
 from cogent3.evolve.models import F81, HKY85, JC69
-from cogent3.evolve.pairwise_distance_numba import (
-    fill_diversity_matrix as numba_fill_diversity_matrix,
-)
 
 warnings.filterwarnings("ignore", "Not using MPI as mpi4py not found")
 
@@ -44,842 +40,96 @@ warnings.filterwarnings("ignore", "Not using MPI as mpi4py not found")
 numpy.seterr(invalid="ignore")
 
 
-class TestPair(TestCase):
-    dna_char_indices = get_moltype_index_array(DNA)
-    rna_char_indices = get_moltype_index_array(RNA)
-    alignment = make_aligned_seqs(
+@pytest.fixture
+def dna_char_indices():
+    """Fixture providing DNA character indices for tests"""
+    return get_moltype_index_array(DNA)
+
+
+@pytest.fixture
+def rna_char_indices():
+    """Fixture providing RNA character indices for tests"""
+    return get_moltype_index_array(RNA)
+
+
+@pytest.fixture
+def basic_alignment():
+    """Fixture providing a basic alignment for tests"""
+    return make_aligned_seqs(
         data=[("s1", "ACGTACGTAC"), ("s2", "GTGTACGTAC")], moltype=DNA
     )
 
-    ambig_alignment = make_aligned_seqs(
+
+@pytest.fixture
+def ambig_alignment():
+    """Fixture providing an ambiguous alignment for tests"""
+    return make_aligned_seqs(
         data=[("s1", "RACGTACGTACN"), ("s2", "AGTGTACGTACA")], moltype=DNA
     )
 
-    diff_alignment = make_aligned_seqs(
+
+@pytest.fixture
+def diff_alignment():
+    """Fixture providing a different alignment for tests"""
+    return make_aligned_seqs(
         data=[("s1", "ACGTACGTTT"), ("s2", "GTGTACGTAC")], moltype=DNA
     )
 
-    def test_char_to_index(self):
-        """should correctly recode a DNA & RNA seqs into indices"""
-        seq = "TCAGRNY?-"
-        expected = [0, 1, 2, 3, -9, -9, -9, -9, -9]
-        indices = seq_to_indices(seq, self.dna_char_indices)
-        assert_equal(indices, expected)
-        seq = "UCAGRNY?-"
-        indices = seq_to_indices(seq, self.rna_char_indices)
-        assert_equal(indices, expected)
-
-    def test_fill_diversity_matrix_all(self):
-        """make correct diversity matrix when all chars valid"""
-        s1 = seq_to_indices("ACGTACGTAC", self.dna_char_indices)
-        s2 = seq_to_indices("GTGTACGTAC", self.dna_char_indices)
-        matrix = numpy.zeros((4, 4), float)
-        # self-self should just be an identity matrix
-        _fill_diversity_matrix(matrix, s1, s1)
-        assert_equal(matrix.sum(), len(s1))
-        assert_equal(
-            matrix,
-            numpy.array(
-                [[2, 0, 0, 0], [0, 3, 0, 0], [0, 0, 3, 0], [0, 0, 0, 2]], float
-            ),
-        )
-
-        # small diffs
-        matrix.fill(0)
-        _fill_diversity_matrix(matrix, s1, s2)
-        assert_equal(
-            matrix,
-            numpy.array(
-                [[2, 0, 0, 0], [1, 2, 0, 0], [0, 0, 2, 1], [0, 0, 0, 2]], float
-            ),
-        )
-
-    def test_fill_diversity_matrix_some(self):
-        """make correct diversity matrix when not all chars valid"""
-        s1 = seq_to_indices("RACGTACGTACN", self.dna_char_indices)
-        s2 = seq_to_indices("AGTGTACGTACA", self.dna_char_indices)
-        matrix = numpy.zeros((4, 4), float)
-        # small diffs
-        matrix.fill(0)
-        _fill_diversity_matrix(matrix, s1, s2)
-        assert_equal(
-            matrix,
-            numpy.array(
-                [[2, 0, 0, 0], [1, 2, 0, 0], [0, 0, 2, 1], [0, 0, 0, 2]], float
-            ),
-        )
-
-    def test_python_vs_numba_fill_matrix(self):
-        """python & cython fill_diversity_matrix give same answer"""
-        s1 = seq_to_indices("RACGTACGTACN", self.dna_char_indices)
-        s2 = seq_to_indices("AGTGTACGTACA", self.dna_char_indices)
-        matrix1 = numpy.zeros((4, 4), float)
-        _fill_diversity_matrix(matrix1, s1, s2)
-        matrix2 = numpy.zeros((4, 4), float)
-        numba_fill_diversity_matrix(matrix2, s1, s2)
-        assert_allclose(matrix1, matrix2)
-
-    def test_hamming_from_matrix(self):
-        """compute hamming from diversity matrix"""
-        s1 = seq_to_indices("ACGTACGTAC", self.dna_char_indices)
-        s2 = seq_to_indices("GTGTACGTAC", self.dna_char_indices)
-        matrix = numpy.zeros((4, 4), float)
-        _fill_diversity_matrix(matrix, s1, s2)
-        total, p, dist, var = _hamming(matrix)
-        self.assertEqual(total, 10.0)
-        self.assertEqual(dist, 2)
-        self.assertEqual(p, 0.2)
-
-    def test_hamming_pair(self):
-        """get distances dict"""
-        calc = HammingPair(DNA, alignment=self.alignment)
-        calc.run(show_progress=False)
-        dists = calc.get_pairwise_distances()
-        dists = dists.to_dict()
-        dist = 2.0
-        expect = {("s1", "s2"): dist, ("s2", "s1"): dist}
-        self.assertEqual(list(dists.keys()), list(expect.keys()))
-        assert_allclose(list(dists.values()), list(expect.values()))
-
-    def test_prop_pair(self):
-        """get distances dict"""
-        calc = ProportionIdenticalPair(DNA, alignment=self.alignment)
-        calc.run(show_progress=False)
-        dists = calc.get_pairwise_distances()
-        dists = dists.to_dict()
-        dist = 0.2
-        expect = {("s1", "s2"): dist, ("s2", "s1"): dist}
-        self.assertEqual(list(dists.keys()), list(expect.keys()))
-        assert_allclose(list(dists.values()), list(expect.values()))
-
-    def test_jc69_from_matrix(self):
-        """compute JC69 from diversity matrix"""
-        s1 = seq_to_indices("ACGTACGTAC", self.dna_char_indices)
-        s2 = seq_to_indices("GTGTACGTAC", self.dna_char_indices)
-        matrix = numpy.zeros((4, 4), float)
-        _fill_diversity_matrix(matrix, s1, s2)
-        total, p, dist, var = _jc69_from_matrix(matrix)
-        self.assertEqual(total, 10.0)
-        self.assertEqual(p, 0.2)
-
-    def test_wrong_moltype(self):
-        """specifying wrong moltype raises ValueError"""
-        with self.assertRaises(ValueError):
-            _ = JC69Pair(PROTEIN, alignment=self.alignment)
-
-    def test_jc69_from_alignment(self):
-        """compute JC69 dists from an alignment"""
-        calc = JC69Pair(DNA, alignment=self.alignment)
-        calc.run(show_progress=False)
-        self.assertEqual(calc.lengths["s1", "s2"], 10)
-        self.assertEqual(calc.proportions["s1", "s2"], 0.2)
-        # value from OSX MEGA 5
-        assert_allclose(calc.dists["s1", "s2"], 0.2326161962)
-        # value**2 from OSX MEGA 5
-        assert_allclose(calc.variances["s1", "s2"], 0.029752066125078681)
-        # value from OSX MEGA 5
-        assert_allclose(calc.stderr["s1", "s2"], 0.1724878724)
-
-        # same answer when using ambiguous alignment
-        calc.run(self.ambig_alignment, show_progress=False)
-        assert_allclose(calc.dists["s1", "s2"], 0.2326161962)
-
-        # but different answer if subsequent alignment is different
-        calc.run(self.diff_alignment, show_progress=False)
-        self.assertTrue(calc.dists["s1", "s2"] != 0.2326161962)
-
-    def test_tn93_from_matrix(self):
-        """compute TN93 distances"""
-        calc = TN93Pair(DNA, alignment=self.alignment)
-        calc.run(show_progress=False)
-        self.assertEqual(calc.lengths["s1", "s2"], 10)
-        self.assertEqual(calc.proportions["s1", "s2"], 0.2)
-        # value from OSX MEGA 5
-        assert_allclose(calc.dists["s1", "s2"], 0.2554128119)
-        # value**2 from OSX MEGA 5
-        assert_allclose(calc.variances["s1", "s2"], 0.04444444445376601)
-        # value from OSX MEGA 5
-        assert_allclose(calc.stderr["s1", "s2"], 0.2108185107)
-
-        # same answer when using ambiguous alignment
-        calc.run(self.ambig_alignment, show_progress=False)
-        assert_allclose(calc.dists["s1", "s2"], 0.2554128119)
-
-        # but different answer if subsequent alignment is different
-        calc.run(self.diff_alignment, show_progress=False)
-        self.assertTrue(calc.dists["s1", "s2"] != 0.2554128119)
-
-    def test_distance_pair(self):
-        """get distances dict"""
-        calc = TN93Pair(DNA, alignment=self.alignment)
-        calc.run(show_progress=False)
-        dists = calc.get_pairwise_distances()
-        dists = dists.to_dict()
-        dist = 0.2554128119
-        expect = {("s1", "s2"): dist, ("s2", "s1"): dist}
-        self.assertEqual(list(dists.keys()), list(expect.keys()))
-        assert_allclose(list(dists.values()), list(expect.values()))
-
-    def test_logdet_pair_dna(self):
-        """logdet should produce distances that match MEGA"""
-        aln = load_aligned_seqs("data/brca1_5.paml", moltype=DNA)
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-        dists = logdet_calc.get_pairwise_distances().to_dict()
-        all_expected = {
-            ("Human", "NineBande"): 0.075336929999999996,
-            ("NineBande", "DogFaced"): 0.0898575452,
-            ("DogFaced", "Human"): 0.1061747919,
-            ("HowlerMon", "DogFaced"): 0.0934480008,
-            ("Mouse", "HowlerMon"): 0.26422862920000001,
-            ("NineBande", "Human"): 0.075336929999999996,
-            ("HowlerMon", "NineBande"): 0.062202897899999998,
-            ("DogFaced", "NineBande"): 0.0898575452,
-            ("DogFaced", "HowlerMon"): 0.0934480008,
-            ("Human", "DogFaced"): 0.1061747919,
-            ("Mouse", "Human"): 0.26539976700000001,
-            ("NineBande", "HowlerMon"): 0.062202897899999998,
-            ("HowlerMon", "Human"): 0.036571181899999999,
-            ("DogFaced", "Mouse"): 0.2652555144,
-            ("HowlerMon", "Mouse"): 0.26422862920000001,
-            ("Mouse", "DogFaced"): 0.2652555144,
-            ("NineBande", "Mouse"): 0.22754789210000001,
-            ("Mouse", "NineBande"): 0.22754789210000001,
-            ("Human", "Mouse"): 0.26539976700000001,
-            ("Human", "HowlerMon"): 0.036571181899999999,
-        }
-        for pair in dists:
-            got = dists[pair]
-            expected = all_expected[pair]
-            assert_allclose(got, expected)
-
-    def test_slice_dmatrix(self):
-        data = {
-            ("ABAYE2984", "Atu3667"): 0.25,
-            ("ABAYE2984", "Avin_42730"): 0.638,
-            ("ABAYE2984", "BAA10469"): None,
-            ("Atu3667", "ABAYE2984"): 0.25,
-            ("Atu3667", "Avin_42730"): 2.368,
-            ("Atu3667", "BAA10469"): 0.25,
-            ("Avin_42730", "ABAYE2984"): 0.638,
-            ("Avin_42730", "Atu3667"): 2.368,
-            ("Avin_42730", "BAA10469"): 1.85,
-            ("BAA10469", "ABAYE2984"): 0.25,
-            ("BAA10469", "Atu3667"): 0.25,
-            ("BAA10469", "Avin_42730"): 1.85,
-        }
-        darr = DistanceMatrix(data)
-        names = darr.template.names[0][:3]
-        got = darr[:3, :3]
-        self.assertEqual(list(got.template.names[0]), names)
-
-    def test_logdet_tk_adjustment(self):
-        """logdet using tamura kumar differs from classic"""
-        aln = load_aligned_seqs("data/brca1_5.paml", moltype=DNA)
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-        tk = logdet_calc.get_pairwise_distances()
-        logdet_calc.run(use_tk_adjustment=False, show_progress=False)
-        not_tk = logdet_calc.get_pairwise_distances()
-        self.assertNotEqual(tk, not_tk)
-
-    def test_logdet_pair_aa(self):
-        """logdet shouldn't fail to produce distances for aa seqs"""
-        aln = load_aligned_seqs("data/brca1_5.paml", moltype=DNA)
-        aln = aln.get_translation()
-        logdet_calc = LogDetPair(moltype=PROTEIN, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-        logdet_calc.get_pairwise_distances()
-
-    def test_logdet_missing_states(self):
-        """should calculate logdet measurement with missing states"""
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTNTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
 
-        dists = logdet_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(list(dists.values())[0] is not None)
-
-        logdet_calc.run(use_tk_adjustment=False, show_progress=False)
-        dists = logdet_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(list(dists.values())[0] is not None)
-
-    def test_logdet_variance(self):
-        """calculate logdet variance consistent with hand calculation"""
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-        self.assertEqual(logdet_calc.variances[1, 1], None)
-
-        index = dict(list(zip("ACGT", list(range(4)))))
-        J = numpy.zeros((4, 4))
-        for p in zip(data[0][1], data[1][1]):
-            J[index[p[0]], index[p[1]]] += 1
-        for i in range(4):
-            if J[i, i] == 0:
-                J[i, i] += 0.5
-        J /= J.sum()
-        M = numpy.linalg.inv(J)
-        var = 0.0
-        for i in range(4):
-            for j in range(4):
-                var += M[j, i] ** 2 * J[i, j] - 1
-        var /= 16 * len(data[0][1])
-
-        logdet_calc.run(use_tk_adjustment=False, show_progress=False)
-        logdet_calc.get_pairwise_distances()
-        assert_allclose(logdet_calc.variances[1, 1], var, atol=1e-3)
-
-    def test_logdet_for_determinant_lte_zero(self):
-        """returns distance of None if the determinant is <= 0"""
-        data = dict(
-            seq1="AGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            seq2="TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-        )
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-        dists = logdet_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(numpy.isnan(list(dists.values())[0]))
-        logdet_calc.run(use_tk_adjustment=False, show_progress=False)
-        dists = logdet_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(numpy.isnan(list(dists.values())[0]))
-
-        # but raises ArithmeticError if told to
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln, invalid_raises=True)
-        with self.assertRaises(ArithmeticError):
-            logdet_calc.run(use_tk_adjustment=True, show_progress=False)
-
-    def test_paralinear_pair_aa(self):
-        """paralinear shouldn't fail to produce distances for aa seqs"""
-        aln = load_aligned_seqs("data/brca1_5.paml", moltype=DNA)
-        aln = aln.get_translation()
-        paralinear_calc = ParalinearPair(moltype=PROTEIN, alignment=aln)
-        paralinear_calc.run(show_progress=False)
-        paralinear_calc.get_pairwise_distances()
-
-    def test_paralinear_distance(self):
-        """calculate paralinear variance consistent with hand calculation"""
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-        paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
-        paralinear_calc.run(show_progress=False)
-
-        index = dict(list(zip("ACGT", list(range(4)))))
-        J = numpy.zeros((4, 4))
-        for p in zip(data[0][1], data[1][1]):
-            J[index[p[0]], index[p[1]]] += 1
-        for i in range(4):
-            if J[i, i] == 0:
-                J[i, i] += 0.5
-        J /= J.sum()
-        numpy.linalg.inv(J)
-        f = J.sum(1), J.sum(0)
-        dist = -0.25 * numpy.log(
-            numpy.linalg.det(J) / numpy.sqrt(f[0].prod() * f[1].prod())
-        )
-
-        assert_allclose(paralinear_calc.dists["seq1", "seq2"], dist)
-
-    def test_paralinear_variance(self):
-        """calculate paralinear variance consistent with hand calculation"""
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-        paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
-        paralinear_calc.run(show_progress=False)
-
-        index = dict(list(zip("ACGT", list(range(4)))))
-        J = numpy.zeros((4, 4))
-        for p in zip(data[0][1], data[1][1]):
-            J[index[p[0]], index[p[1]]] += 1
-        for i in range(4):
-            if J[i, i] == 0:
-                J[i, i] += 0.5
-        J /= J.sum()
-        M = numpy.linalg.inv(J)
-        f = J.sum(1), J.sum(0)
-        var = 0.0
-        for i in range(4):
-            for j in range(4):
-                var += M[j, i] ** 2 * J[i, j]
-            var -= 1 / numpy.sqrt(f[0][i] * f[1][i])
-        var /= 16 * len(data[0][1])
-
-        assert_allclose(paralinear_calc.variances[1, 1], var, atol=1e-3)
-
-    def test_paralinear_for_determinant_lte_zero(self):
-        """returns distance of None if the determinant is <= 0"""
-        data = dict(
-            seq1="AGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            seq2="TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-        )
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-
-        paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
-        paralinear_calc.run(show_progress=False)
-        dists = paralinear_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(numpy.isnan(list(dists.values())[0]))
-        paralinear_calc.run(show_progress=False)
-        dists = paralinear_calc.get_pairwise_distances().to_dict()
-        self.assertTrue(numpy.isnan(list(dists.values())[0]))
-
-    def test_paralinear_pair_dna(self):
-        """calculate paralinear distance consistent with logdet distance"""
-        data = [
-            (
-                "seq1",
-                "TAATTCATTGGGACGTCGAATCCGGCAGTCCTGCCGCAAAAGCTTCCGGAATCGAATTTTGGCA",
-            ),
-            (
-                "seq2",
-                "AAAAAAAAAAAAAAAACCCCCCCCCCCCCCCCTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGG",
-            ),
-        ]
-        aln = make_aligned_seqs(data=data, moltype=DNA)
-        paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
-        paralinear_calc.run(show_progress=False)
-        logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
-        logdet_calc.run(show_progress=False)
-        self.assertEqual(logdet_calc.dists[1, 1], paralinear_calc.dists[1, 1])
-        self.assertEqual(paralinear_calc.variances[1, 1], logdet_calc.variances[1, 1])
-
-    def test_duplicated(self):
-        """correctly identifies duplicates"""
-
-        def get_calc(data):
-            aln = make_aligned_seqs(data=data, moltype=DNA)
-            calc = ParalinearPair(moltype=DNA, alignment=aln)
-            calc(show_progress=False)
-            return calc
-
-        # no duplicates
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        calc = get_calc(data)
-        self.assertEqual(calc.duplicated, None)
-        data = [
-            (
-                "seq1",
-                "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
-            ),
-            (
-                "seq2",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-            (
-                "seq3",
-                "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
-            ),
-        ]
-        calc = get_calc(data)
-        self.assertTrue(
-            {"seq2": ["seq3"]} == calc.duplicated
-            or {"seq3": ["seq2"]} == calc.duplicated
-        )
-        # default to get all pairwise distances
-        pwds = calc.get_pairwise_distances().to_dict()
-        self.assertEqual(pwds[("seq2", "seq3")], 0.0)
-        self.assertEqual(pwds[("seq2", "seq1")], pwds[("seq3", "seq1")])
-
-        # only unique seqs when using include_duplicates=False
-
-        pwds = calc.get_pairwise_distances(include_duplicates=False).to_dict()
-        present = list(calc.duplicated.keys())[0]
-        missing = calc.duplicated[present][0]
-        self.assertEqual(set([(present, missing)]), set([("seq2", "seq3")]))
-        self.assertTrue((present, "seq1") in pwds)
-        self.assertFalse((missing, "seq1") in pwds)
-
-
-class TestGetDisplayCalculators(TestCase):
-    def test_get_calculator(self):
-        """exercising getting specified calculator"""
-        for key in _calculators:
-            get_distance_calculator(key)
-            get_distance_calculator(key.upper())
-
-        with self.assertRaises(ValueError):
-            get_distance_calculator("blahblah")
-
-    def test_available_distances(self):
-        """available_distances has correct content"""
-        content = available_distances()
-        self.assertEqual(content.shape, (6, 2))
-        self.assertEqual(content["tn93", 1], "dna, rna")
-
-
-class TestDistanceMatrix(TestCase):
-    def test_to_dict(self):
-        """distance matrix correctly produces a 1D dict"""
-        data = {("s1", "s2"): 0.25, ("s2", "s1"): 0.25}
-        dmat = DistanceMatrix(data)
-        got = dmat.to_dict()
-        self.assertEqual(got, data)
-
-    def test_matrix_dtype(self):
-        """tests DistanceMatrix correctly accepts the data with proper dtype"""
-        data = {
-            ("ABAYE2984", "Atu3667"): None,
-            ("ABAYE2984", "Avin_42730"): 0.638,
-            ("ABAYE2984", "BAA10469"): None,
-            ("Atu3667", "ABAYE2984"): None,
-            ("Atu3667", "Avin_42730"): 2.368,
-            ("Atu3667", "BAA10469"): None,
-            ("Avin_42730", "ABAYE2984"): 0.638,
-            ("Avin_42730", "Atu3667"): 2.368,
-            ("Avin_42730", "BAA10469"): 1.85,
-            ("BAA10469", "ABAYE2984"): None,
-            ("BAA10469", "Atu3667"): None,
-            ("BAA10469", "Avin_42730"): 1.85,
-        }
-        names = set()
-        for p in data:
-            names.update(p)
-
-        # tests when data has None values and DistanceMatrix using dtype('float')
-        darr = DistanceMatrix(data)
-        self.assertEqual(darr.shape, (4, 4))
-        self.assertEqual(set(darr.names), names)
-        for (a, b), dist in data.items():
-            if dist is None:
-                assert numpy.isnan(darr[a, b])
-            else:
-                assert_allclose(dist, darr[a, b])
-
-        data = {
-            ("ABAYE2984", "Atu3667"): "None",
-            ("ABAYE2984", "Avin_42730"): 0.638,
-            ("ABAYE2984", "BAA10469"): None,
-            ("Atu3667", "ABAYE2984"): None,
-            ("Atu3667", "Avin_42730"): 2.368,
-            ("Atu3667", "BAA10469"): "None",
-            ("Avin_42730", "ABAYE2984"): 0.638,
-            ("Avin_42730", "Atu3667"): 2.368,
-            ("Avin_42730", "BAA10469"): 1.85,
-            ("BAA10469", "ABAYE2984"): None,
-            ("BAA10469", "Atu3667"): None,
-            ("BAA10469", "Avin_42730"): 1.85,
-        }
+@pytest.fixture
+def alignment():
+    return make_aligned_seqs(
+        data=[("s1", "ACGTACGTAC"), ("s2", "GTGTACGTAC")], moltype=DNA
+    )
 
-        # tests when data has str values and DistanceMatrix using dtype('float')
-        with self.assertRaises(ValueError):
-            darr = DistanceMatrix(data)
-
-    def test_take_dists(self):
-        """subsets the distance matrix"""
-        data = {
-            ("ABAYE2984", "Atu3667"): 0.25,
-            ("ABAYE2984", "Avin_42730"): 0.638,
-            ("ABAYE2984", "BAA10469"): None,
-            ("Atu3667", "ABAYE2984"): 0.25,
-            ("Atu3667", "Avin_42730"): 2.368,
-            ("Atu3667", "BAA10469"): 0.25,
-            ("Avin_42730", "ABAYE2984"): 0.638,
-            ("Avin_42730", "Atu3667"): 2.368,
-            ("Avin_42730", "BAA10469"): 1.85,
-            ("BAA10469", "ABAYE2984"): 0.25,
-            ("BAA10469", "Atu3667"): 0.25,
-            ("BAA10469", "Avin_42730"): 1.85,
-        }
-        darr = DistanceMatrix(data)
-        got1 = darr.take_dists(["ABAYE2984", "Atu3667", "Avin_42730"])
-        got2 = darr.take_dists("BAA10469", negate=True)
-        assert_allclose(got1.array.astype(float), got2.array.astype(float))
-
-    def test_build_phylogeny(self):
-        """build a NJ tree"""
-        from cogent3 import make_tree
-
-        dists = {
-            ("DogFaced", "FlyingFox"): 0.05,
-            ("DogFaced", "FreeTaile"): 0.14,
-            ("DogFaced", "LittleBro"): 0.16,
-            ("DogFaced", "TombBat"): 0.15,
-            ("FlyingFox", "DogFaced"): 0.05,
-            ("FlyingFox", "FreeTaile"): 0.12,
-            ("FlyingFox", "LittleBro"): 0.13,
-            ("FlyingFox", "TombBat"): 0.14,
-            ("FreeTaile", "DogFaced"): 0.14,
-            ("FreeTaile", "FlyingFox"): 0.12,
-            ("FreeTaile", "LittleBro"): 0.09,
-            ("FreeTaile", "TombBat"): 0.1,
-            ("LittleBro", "DogFaced"): 0.16,
-            ("LittleBro", "FlyingFox"): 0.13,
-            ("LittleBro", "FreeTaile"): 0.09,
-            ("LittleBro", "TombBat"): 0.12,
-            ("TombBat", "DogFaced"): 0.15,
-            ("TombBat", "FlyingFox"): 0.14,
-            ("TombBat", "FreeTaile"): 0.1,
-            ("TombBat", "LittleBro"): 0.12,
-        }
-        dists = DistanceMatrix(dists)
-        got = dists.quick_tree(show_progress=False)
-        expect = make_tree(
-            treestring="((TombBat,(DogFaced,FlyingFox)),LittleBro,FreeTaile)"
-        )
-        self.assertTrue(expect.same_topology(got))
-
-    def test_names(self):
-        """names property works"""
-        data = {
-            ("ABAYE2984", "Atu3667"): 0.25,
-            ("ABAYE2984", "Avin_42730"): 0.638,
-            ("ABAYE2984", "BAA10469"): None,
-            ("Atu3667", "ABAYE2984"): 0.25,
-            ("Atu3667", "Avin_42730"): 2.368,
-            ("Atu3667", "BAA10469"): 0.25,
-            ("Avin_42730", "ABAYE2984"): 0.638,
-            ("Avin_42730", "Atu3667"): 2.368,
-            ("Avin_42730", "BAA10469"): 1.85,
-            ("BAA10469", "ABAYE2984"): 0.25,
-            ("BAA10469", "Atu3667"): 0.25,
-            ("BAA10469", "Avin_42730"): 1.85,
-        }
-        names = set()
-        for p in data:
-            names.update(p)
-        darr = DistanceMatrix(data)
-        self.assertEqual(set(darr.names), names)
-        darr = darr.drop_invalid()
-        for n in ("ABAYE2984", "BAA10469"):
-            names.remove(n)
-        self.assertEqual(set(darr.names), names)
-
-
-class DistancesTests(TestCase):
-    def setUp(self):
-        self.al = make_aligned_seqs(
-            data={
-                "a": "GTACGTACGATC",
-                "b": "GTACGTACGTAC",
-                "c": "GTACGTACGTTC",
-                "e": "GTACGTACTGGT",
-            },
-            moltype="dna",
-        )
-        self.collection = make_unaligned_seqs(
-            data={
-                "a": "GTACGTACGATC",
-                "b": "GTACGTACGTAC",
-                "c": "GTACGTACGTTC",
-                "e": "GTACGTACTGGT",
-            },
-            moltype="dna",
-        )
-
-    def assertDistsAlmostEqual(self, expected, observed, precision=4):
-        observed = dict([(frozenset(k), v) for (k, v) in list(observed.items())])
-        expected = dict([(frozenset(k), v) for (k, v) in list(expected.items())])
-        for key in expected:
-            self.assertAlmostEqual(expected[key], observed[key], precision)
-
-    def test_EstimateDistances(self):
-        """testing (well, exercising at least), EstimateDistances"""
-        d = EstimateDistances(self.al, JC69())
-        d.run(show_progress=False)
-        canned_result = {
-            ("b", "e"): 0.440840,
-            ("c", "e"): 0.440840,
-            ("a", "c"): 0.088337,
-            ("a", "b"): 0.188486,
-            ("a", "e"): 0.440840,
-            ("b", "c"): 0.0883373,
-        }
-        result = d.get_pairwise_distances().to_dict()
-        self.assertDistsAlmostEqual(canned_result, result)
-
-        # excercise writing to file
-        d.write("junk.txt")
-        try:
-            os.remove("junk.txt")
-        except OSError:
-            pass  # probably parallel
-
-    def test_EstimateDistancesWithMotifProbs(self):
-        """EstimateDistances with supplied motif probs"""
-        motif_probs = {"A": 0.1, "C": 0.2, "G": 0.2, "T": 0.5}
-        d = EstimateDistances(self.al, HKY85(), motif_probs=motif_probs)
-        d.run(show_progress=False)
-        canned_result = {
-            ("a", "c"): 0.07537,
-            ("b", "c"): 0.07537,
-            ("a", "e"): 0.39921,
-            ("a", "b"): 0.15096,
-            ("b", "e"): 0.39921,
-            ("c", "e"): 0.37243,
-        }
-        result = d.get_pairwise_distances().to_dict()
-        self.assertDistsAlmostEqual(canned_result, result)
-
-    def test_EstimateDistances_fromThreeway(self):
-        """testing (well, exercising at least), EsimateDistances fromThreeway"""
-        d = EstimateDistances(self.al, JC69(), threeway=True)
-        d.run(show_progress=False)
-        canned_result = {
-            ("b", "e"): 0.495312,
-            ("c", "e"): 0.479380,
-            ("a", "c"): 0.089934,
-            ("a", "b"): 0.190021,
-            ("a", "e"): 0.495305,
-            ("b", "c"): 0.0899339,
-        }
-        result = d.get_pairwise_distances(summary_function="mean").to_dict()
-        self.assertDistsAlmostEqual(canned_result, result)
-
-    def test_EstimateDistances_fromUnaligned(self):
-        """Excercising estimate distances from unaligned sequences"""
-        d = EstimateDistances(
-            self.collection, JC69(), do_pair_align=True, rigorous_align=True
-        )
-        d.run(show_progress=False)
-        canned_result = {
-            ("b", "e"): 0.440840,
-            ("c", "e"): 0.440840,
-            ("a", "c"): 0.088337,
-            ("a", "b"): 0.188486,
-            ("a", "e"): 0.440840,
-            ("b", "c"): 0.0883373,
-        }
-        result = d.get_pairwise_distances().to_dict()
-        self.assertDistsAlmostEqual(canned_result, result)
-
-        d = EstimateDistances(
-            self.collection, JC69(), do_pair_align=True, rigorous_align=False
-        )
-        d.run(show_progress=False)
-        canned_result = {
-            ("b", "e"): 0.440840,
-            ("c", "e"): 0.440840,
-            ("a", "c"): 0.088337,
-            ("a", "b"): 0.188486,
-            ("a", "e"): 0.440840,
-            ("b", "c"): 0.0883373,
-        }
-        result = d.get_pairwise_distances().to_dict()
-        self.assertDistsAlmostEqual(canned_result, result)
-
-    def test_EstimateDistances_other_model_params(self):
-        """test getting other model params from EstimateDistances"""
-        d = EstimateDistances(self.al, HKY85(), est_params=["kappa"])
-        d.run(show_progress=False)
-        # this will be a Number object with Mean, Median etc ..
-        kappa = d.get_param_values("kappa")
-        self.assertAlmostEqual(kappa.mean, 0.8939, 4)
-        # this will be a dict with pairwise instances, it's called by the above
-        # method, so the correctness of it's values is already checked
-        kappa = d.get_pairwise_param("kappa")
-
-    def test_EstimateDistances_modify_lf(self):
-        """tests modifying the lf"""
-
-        def constrain_fit(lf):
-            lf.set_param_rule("kappa", is_constant=True)
-            lf.optimise(local=True, show_progress=False)
-            return lf
-
-        d = EstimateDistances(self.al, HKY85(), modify_lf=constrain_fit)
-        d.run(show_progress=False)
-        result = d.get_pairwise_distances().to_dict()
-        d = EstimateDistances(self.al, F81())
-        d.run(show_progress=False)
-        expect = d.get_pairwise_distances().to_dict()
-        self.assertDistsAlmostEqual(expect, result)
-
-    def test_get_raw_estimates(self):
-        """correctly return raw result object"""
-        d = EstimateDistances(self.al, HKY85(), est_params=["kappa"])
-        d.run(show_progress=False)
-        expect = {
-            ("a", "b"): {
-                "kappa": 1.0000226766004808e-06,
-                "length": 0.18232155856115662,
-            },
-            ("a", "c"): {
-                "kappa": 1.0010380037049357e-06,
-                "length": 0.087070406623635604,
-            },
-            ("a", "e"): {"kappa": 2.3965871843412687, "length": 0.4389176272584539},
-            ("b", "e"): {"kappa": 2.3965871854366592, "length": 0.43891762729173389},
-            ("b", "c"): {
-                "kappa": 1.0010380037049357e-06,
-                "length": 0.087070406623635604,
-            },
-            ("c", "e"): {"kappa": 0.57046787478038707, "length": 0.43260232210282784},
-        }
-        got = d.get_all_param_values()
-        for pair in expect:
-            for param in expect[pair]:
-                self.assertAlmostEqual(got[pair][param], expect[pair][param], places=6)
-
-    def test_no_calc(self):
-        """returns None if no calculation done"""
-        al = load_aligned_seqs("data/brca1_5.paml")
-        d = EstimateDistances(al, submodel=HKY85())
-        self.assertEqual(d.get_pairwise_distances(), None)
-
-    def test_to_table(self):
-        """converts a distance matrix to a Table"""
-        data = {
-            ("A", "B"): 2,
-            ("A", "C"): 3,
-            ("B", "C"): 1,
-            ("B", "A"): 2,
-            ("C", "A"): 3,
-            ("C", "B"): 1,
-        }
-        darr = DistanceMatrix(data)
-        table = darr.to_table()
-        self.assertEqual(table.shape, (3, 4))
-        self.assertEqual(table.columns["names"].tolist(), list(darr.names))
-        self.assertEqual(table["A", "B"], 2)
-        self.assertEqual(table["A", "A"], 0)
+
+@pytest.fixture(scope="session")
+def brca1_5(DATA_DIR):
+    return load_aligned_seqs(DATA_DIR / "brca1_5.paml", moltype=DNA)
+
+
+def test_char_to_index(dna_char_indices, rna_char_indices):
+    """should correctly recode a DNA & RNA seqs into indices"""
+    seq = "TCAGRNY?-"
+    expected = [0, 1, 2, 3, -9, -9, -9, -9, -9]
+    indices = seq_to_indices(seq, dna_char_indices)
+    assert_equal(indices, expected)
+    seq = "UCAGRNY?-"
+    indices = seq_to_indices(seq, rna_char_indices)
+    assert_equal(indices, expected)
+
+
+@pytest.fixture(scope="session")
+def al():
+    return make_aligned_seqs(
+        data={
+            "a": "GTACGTACGATC",
+            "b": "GTACGTACGTAC",
+            "c": "GTACGTACGTTC",
+            "e": "GTACGTACTGGT",
+        },
+        moltype="dna",
+    )
+
+
+@pytest.fixture(scope="session")
+def collection():
+    return make_unaligned_seqs(
+        data={
+            "a": "GTACGTACGATC",
+            "b": "GTACGTACGTAC",
+            "c": "GTACGTACGTTC",
+            "e": "GTACGTACTGGT",
+        },
+        moltype="dna",
+    )
+
+
+def assert_dists_approx_equal(expected, observed, atol=1e-6):
+    observed = dict([(frozenset(k), v) for (k, v) in list(observed.items())])
+    expected = dict([(frozenset(k), v) for (k, v) in list(expected.items())])
+    for key in expected:
+        assert_allclose(expected[key], observed[key], rtol=1e-6, atol=atol)
 
 
 @pytest.fixture
@@ -1010,3 +260,930 @@ def test_distance_matrix_from_array_names():
     names = list("abcd")
     got = DistanceMatrix.from_array_names(data, names)
     assert got["b", "c"] == 0.1
+
+
+def test_fill_diversity_matrix_all(dna_char_indices):
+    """make correct diversity matrix when all chars valid"""
+    s1 = seq_to_indices("ACGTACGTAC", dna_char_indices)
+    s2 = seq_to_indices("GTGTACGTAC", dna_char_indices)
+    matrix = numpy.zeros((4, 4), float)
+    # self-self should just be an identity matrix
+    fill_diversity_matrix(matrix, s1, s1)
+    assert_equal(matrix.sum(), len(s1))
+    assert_equal(
+        matrix,
+        numpy.array([[2, 0, 0, 0], [0, 3, 0, 0], [0, 0, 3, 0], [0, 0, 0, 2]], float),
+    )
+
+    # small diffs
+    matrix.fill(0)
+    fill_diversity_matrix(matrix, s1, s2)
+    assert_equal(
+        matrix,
+        numpy.array([[2, 0, 0, 0], [1, 2, 0, 0], [0, 0, 2, 1], [0, 0, 0, 2]], float),
+    )
+
+
+def test_fill_diversity_matrix_some(dna_char_indices):
+    """make correct diversity matrix when not all chars valid"""
+    s1 = seq_to_indices("RACGTACGTACN", dna_char_indices)
+    s2 = seq_to_indices("AGTGTACGTACA", dna_char_indices)
+    matrix = numpy.zeros((4, 4), float)
+    # self-self should just be an identity matrix
+    fill_diversity_matrix(matrix, s1, s1)
+    assert_equal(matrix.sum(), 10)  # 2 invalid chars
+    assert_equal(
+        matrix,
+        numpy.array([[2, 0, 0, 0], [0, 3, 0, 0], [0, 0, 3, 0], [0, 0, 0, 2]], float),
+    )
+
+    # small diffs
+    matrix.fill(0)
+    fill_diversity_matrix(matrix, s1, s2)
+    assert_equal(
+        matrix,
+        numpy.array([[2, 0, 0, 0], [1, 2, 0, 0], [0, 0, 2, 1], [0, 0, 0, 2]], float),
+    )
+
+
+def test_hamming_from_matrix(dna_char_indices):
+    """compute hamming from diversity matrix"""
+    s1 = seq_to_indices("ACGTACGTAC", dna_char_indices)
+    s2 = seq_to_indices("GTGTACGTAC", dna_char_indices)
+    matrix = numpy.zeros((4, 4), float)
+    fill_diversity_matrix(matrix, s1, s2)
+    total = 0
+    for i in range(4):
+        total += matrix[i, i]
+    assert_equal(total / len(s1), 0.8)
+
+
+def test_hamming_pair(basic_alignment):
+    """get distances dict"""
+    calc = HammingPair(DNA, alignment=basic_alignment)
+    calc.run(show_progress=False)
+    dists = calc.get_pairwise_distances()
+    assert_equal(dists.to_dict()["s1", "s2"], 2)
+
+
+def test_prop_pair(basic_alignment):
+    """get distances dict"""
+    calc = ProportionIdenticalPair(DNA, alignment=basic_alignment)
+    calc.run(show_progress=False)
+    dists = calc.get_pairwise_distances()
+    assert_equal(dists.to_dict()["s1", "s2"], 0.2)
+
+
+def test_jc69_from_matrix(dna_char_indices):
+    """compute JC69 from diversity matrix"""
+    s1 = seq_to_indices("ACGTACGTAC", dna_char_indices)
+    s2 = seq_to_indices("GTGTACGTAC", dna_char_indices)
+    matrix = numpy.zeros((4, 4), float)
+    fill_diversity_matrix(matrix, s1, s2)
+    total = 0
+    for i in range(4):
+        total += matrix[i, i]
+    p = 1 - (total / len(s1))
+    assert_allclose(p, 0.2)
+
+
+def test_wrong_moltype(basic_alignment):
+    """specifying wrong moltype raises ValueError"""
+    with pytest.raises(ValueError):
+        _ = JC69Pair(PROTEIN, alignment=basic_alignment)
+
+
+def test_jc69_from_alignment(basic_alignment, ambig_alignment, diff_alignment):
+    """compute JC69 dists from an alignment"""
+    calc = JC69Pair(DNA, alignment=basic_alignment)
+    calc.run(show_progress=False)
+    assert_equal(calc.lengths["s1", "s2"], 10)
+    assert_equal(calc.proportions["s1", "s2"], 0.2)
+    # value from OSX MEGA 5
+    assert_allclose(calc.dists["s1", "s2"], 0.2326161962)
+    # value**2 from OSX MEGA 5
+    assert_allclose(calc.variances["s1", "s2"], 0.029752066125078681)
+    # value from OSX MEGA 5
+    assert_allclose(calc.stderr["s1", "s2"], 0.1724878724)
+
+    # same answer when using ambiguous alignment
+    calc.run(ambig_alignment, show_progress=False)
+    assert_allclose(calc.dists["s1", "s2"], 0.2326161962)
+
+    # but different answer if subsequent alignment is different
+    calc.run(diff_alignment, show_progress=False)
+    assert calc.dists["s1", "s2"] != 0.2326161962
+
+
+def test_tn93_from_matrix(basic_alignment, ambig_alignment, diff_alignment):
+    """compute TN93 distances"""
+    calc = TN93Pair(DNA, alignment=basic_alignment)
+    calc.run(show_progress=False)
+    assert_equal(calc.lengths["s1", "s2"], 10)
+    assert_equal(calc.proportions["s1", "s2"], 0.2)
+    # value from OSX MEGA 5
+    assert_allclose(calc.dists["s1", "s2"], 0.2554128119)
+    # value**2 from OSX MEGA 5
+    assert_allclose(calc.variances["s1", "s2"], 0.04444444445376601)
+    # value from OSX MEGA 5
+    assert_allclose(calc.stderr["s1", "s2"], 0.2108185107)
+
+    # same answer when using ambiguous alignment
+    calc.run(ambig_alignment, show_progress=False)
+    assert_allclose(calc.dists["s1", "s2"], 0.2554128119)
+
+    # but different answer if subsequent alignment is different
+    calc.run(diff_alignment, show_progress=False)
+    assert calc.dists["s1", "s2"] != 0.2554128119
+
+
+def test_distance_pair(alignment):
+    """get distances dict"""
+    calc = TN93Pair(DNA, alignment=alignment)
+    calc.run(show_progress=False)
+    dists = calc.get_pairwise_distances()
+    dists = dists.to_dict()
+    dist = 0.2554128119
+    expect = {("s1", "s2"): dist, ("s2", "s1"): dist}
+    assert list(dists.keys()) == list(expect.keys())
+    assert_allclose(list(dists.values()), list(expect.values()))
+
+
+def test_logdet_pair_dna(brca1_5):
+    """logdet should produce distances that match MEGA"""
+    aln = brca1_5
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+    dists = logdet_calc.get_pairwise_distances().to_dict()
+    all_expected = {
+        ("Human", "NineBande"): 0.075336929999999996,
+        ("NineBande", "DogFaced"): 0.0898575452,
+        ("DogFaced", "Human"): 0.1061747919,
+        ("HowlerMon", "DogFaced"): 0.0934480008,
+        ("Mouse", "HowlerMon"): 0.26422862920000001,
+        ("NineBande", "Human"): 0.075336929999999996,
+        ("HowlerMon", "NineBande"): 0.062202897899999998,
+        ("DogFaced", "NineBande"): 0.0898575452,
+        ("DogFaced", "HowlerMon"): 0.0934480008,
+        ("Human", "DogFaced"): 0.1061747919,
+        ("Mouse", "Human"): 0.26539976700000001,
+        ("NineBande", "HowlerMon"): 0.062202897899999998,
+        ("HowlerMon", "Human"): 0.036571181899999999,
+        ("DogFaced", "Mouse"): 0.2652555144,
+        ("HowlerMon", "Mouse"): 0.26422862920000001,
+        ("Mouse", "DogFaced"): 0.2652555144,
+        ("NineBande", "Mouse"): 0.22754789210000001,
+        ("Mouse", "NineBande"): 0.22754789210000001,
+        ("Human", "Mouse"): 0.26539976700000001,
+        ("Human", "HowlerMon"): 0.036571181899999999,
+    }
+    for pair in dists:
+        got = dists[pair]
+        expected = all_expected[pair]
+        assert_allclose(got, expected)
+
+
+def test_slice_dmatrix():
+    data = {
+        ("ABAYE2984", "Atu3667"): 0.25,
+        ("ABAYE2984", "Avin_42730"): 0.638,
+        ("ABAYE2984", "BAA10469"): None,
+        ("Atu3667", "ABAYE2984"): 0.25,
+        ("Atu3667", "Avin_42730"): 2.368,
+        ("Atu3667", "BAA10469"): 0.25,
+        ("Avin_42730", "ABAYE2984"): 0.638,
+        ("Avin_42730", "Atu3667"): 2.368,
+        ("Avin_42730", "BAA10469"): 1.85,
+        ("BAA10469", "ABAYE2984"): 0.25,
+        ("BAA10469", "Atu3667"): 0.25,
+        ("BAA10469", "Avin_42730"): 1.85,
+    }
+    darr = DistanceMatrix(data)
+    names = darr.template.names[0][:3]
+    got = darr[:3, :3]
+    assert list(got.template.names[0]) == names
+
+
+def test_logdet_tk_adjustment(brca1_5):
+    """logdet using tamura kumar differs from classic"""
+    aln = brca1_5
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+    tk = logdet_calc.get_pairwise_distances()
+    logdet_calc.run(use_tk_adjustment=False, show_progress=False)
+    not_tk = logdet_calc.get_pairwise_distances()
+    assert tk != not_tk
+
+
+def test_logdet_pair_aa(brca1_5):
+    """logdet shouldn't fail to produce distances for aa seqs"""
+    aln = brca1_5
+    aln = aln.get_translation()
+    logdet_calc = LogDetPair(moltype=PROTEIN, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+    dists = logdet_calc.get_pairwise_distances()
+    assert isinstance(dists, DistanceMatrix)
+
+
+def test_logdet_missing_states():
+    """should calculate logdet measurement with missing states"""
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTNTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+
+    dists = logdet_calc.get_pairwise_distances().to_dict()
+    assert list(dists.values())[0] is not None
+
+    logdet_calc.run(use_tk_adjustment=False, show_progress=False)
+    dists = logdet_calc.get_pairwise_distances().to_dict()
+    assert list(dists.values())[0] is not None
+
+
+def test_logdet_variance():
+    """calculate logdet variance consistent with hand calculation"""
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+    assert logdet_calc.variances[1, 1] is None
+
+    index = dict(list(zip("ACGT", list(range(4)))))
+    J = numpy.zeros((4, 4))
+    for p in zip(data[0][1], data[1][1]):
+        J[index[p[0]], index[p[1]]] += 1
+    for i in range(4):
+        if J[i, i] == 0:
+            J[i, i] += 0.5
+    J /= J.sum()
+    M = numpy.linalg.inv(J)
+    var = 0.0
+    for i in range(4):
+        for j in range(4):
+            var += M[j, i] ** 2 * J[i, j] - 1
+    var /= 16 * len(data[0][1])
+
+    logdet_calc.run(use_tk_adjustment=False, show_progress=False)
+    logdet_calc.get_pairwise_distances()
+    assert_allclose(logdet_calc.variances[1, 1], var, atol=1e-3)
+
+
+def test_logdet_for_determinant_lte_zero():
+    """returns distance of None if the determinant is <= 0"""
+    data = dict(
+        seq1="AGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        seq2="TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+    )
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+    dists = logdet_calc.get_pairwise_distances().to_dict()
+    assert numpy.isnan(list(dists.values())[0])
+    logdet_calc.run(use_tk_adjustment=False, show_progress=False)
+    dists = logdet_calc.get_pairwise_distances().to_dict()
+    assert numpy.isnan(list(dists.values())[0])
+
+    # but raises ArithmeticError if told to
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln, invalid_raises=True)
+    with pytest.raises(ArithmeticError):
+        logdet_calc.run(use_tk_adjustment=True, show_progress=False)
+
+
+def test_paralinear_pair_aa():
+    """paralinear shouldn't fail to produce distances for aa seqs"""
+    aln = load_aligned_seqs("data/brca1_5.paml", moltype=DNA)
+    aln = aln.get_translation()
+    paralinear_calc = ParalinearPair(moltype=PROTEIN, alignment=aln)
+    paralinear_calc.run(show_progress=False)
+    dists = paralinear_calc.get_pairwise_distances()
+    assert isinstance(dists, DistanceMatrix)
+
+
+def test_paralinear_distance():
+    """calculate paralinear variance consistent with hand calculation"""
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
+    paralinear_calc.run(show_progress=False)
+
+    index = dict(list(zip("ACGT", list(range(4)))))
+    J = numpy.zeros((4, 4))
+    for p in zip(data[0][1], data[1][1]):
+        J[index[p[0]], index[p[1]]] += 1
+    for i in range(4):
+        if J[i, i] == 0:
+            J[i, i] += 0.5
+    J /= J.sum()
+    numpy.linalg.inv(J)
+    f = J.sum(1), J.sum(0)
+    dist = -0.25 * numpy.log(
+        numpy.linalg.det(J) / numpy.sqrt(f[0].prod() * f[1].prod())
+    )
+
+    assert_allclose(paralinear_calc.dists["seq1", "seq2"], dist)
+
+
+def test_paralinear_variance():
+    """calculate paralinear variance consistent with hand calculation"""
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
+    paralinear_calc.run(show_progress=False)
+
+    index = dict(list(zip("ACGT", list(range(4)))))
+    J = numpy.zeros((4, 4))
+    for p in zip(data[0][1], data[1][1]):
+        J[index[p[0]], index[p[1]]] += 1
+    for i in range(4):
+        if J[i, i] == 0:
+            J[i, i] += 0.5
+    J /= J.sum()
+    M = numpy.linalg.inv(J)
+    f = J.sum(1), J.sum(0)
+    var = 0.0
+    for i in range(4):
+        for j in range(4):
+            var += M[j, i] ** 2 * J[i, j]
+        var -= 1 / numpy.sqrt(f[0][i] * f[1][i])
+    var /= 16 * len(data[0][1])
+
+    assert_allclose(paralinear_calc.variances[1, 1], var, atol=1e-3)
+
+
+def test_paralinear_for_determinant_lte_zero():
+    """returns distance of None if the determinant is <= 0"""
+    data = dict(
+        seq1="AGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        seq2="TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+    )
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+
+    paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
+    paralinear_calc.run(show_progress=False)
+    dists = paralinear_calc.get_pairwise_distances().to_dict()
+    assert numpy.isnan(list(dists.values())[0])
+    paralinear_calc.run(show_progress=False)
+    dists = paralinear_calc.get_pairwise_distances().to_dict()
+    assert numpy.isnan(list(dists.values())[0])
+
+
+def test_paralinear_pair_dna():
+    """calculate paralinear distance consistent with logdet distance"""
+    data = [
+        (
+            "seq1",
+            "TAATTCATTGGGACGTCGAATCCGGCAGTCCTGCCGCAAAAGCTTCCGGAATCGAATTTTGGCA",
+        ),
+        (
+            "seq2",
+            "AAAAAAAAAAAAAAAACCCCCCCCCCCCCCCCTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGG",
+        ),
+    ]
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    paralinear_calc = ParalinearPair(moltype=DNA, alignment=aln)
+    paralinear_calc.run(show_progress=False)
+    logdet_calc = LogDetPair(moltype=DNA, alignment=aln)
+    logdet_calc.run(show_progress=False)
+    assert logdet_calc.dists[1, 1] == paralinear_calc.dists[1, 1]
+    assert paralinear_calc.variances[1, 1] == logdet_calc.variances[1, 1]
+
+
+def get_calc(data):
+    aln = make_aligned_seqs(data=data, moltype=DNA)
+    calc = ParalinearPair(moltype=DNA, alignment=aln)
+    calc(show_progress=False)
+    return calc
+
+
+def test_duplicated():
+    """correctly identifies duplicates"""
+
+    # no duplicates
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    calc = get_calc(data)
+    assert calc.duplicated is None
+    data = [
+        (
+            "seq1",
+            "GGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        ),
+        (
+            "seq2",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+        (
+            "seq3",
+            "TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+        ),
+    ]
+    calc = get_calc(data)
+    assert {"seq2": ["seq3"]} == calc.duplicated or {
+        "seq3": ["seq2"]
+    } == calc.duplicated
+    # default to get all pairwise distances
+    pwds = calc.get_pairwise_distances().to_dict()
+    assert pwds[("seq2", "seq3")] == 0.0
+    assert pwds[("seq2", "seq1")] == pwds[("seq3", "seq1")]
+
+    # only unique seqs when using include_duplicates=False
+
+    pwds = calc.get_pairwise_distances(include_duplicates=False).to_dict()
+    present = list(calc.duplicated.keys())[0]
+    missing = calc.duplicated[present][0]
+    assert {(present, missing)} == {("seq2", "seq3")}
+    assert (present, "seq1") in pwds
+    assert (missing, "seq1") not in pwds
+
+
+def test_get_calculator():
+    """exercising getting specified calculator"""
+    for key in _calculators:
+        get_distance_calculator(key)
+        get_distance_calculator(key.upper())
+
+    with pytest.raises(ValueError):
+        get_distance_calculator("blahblah")
+
+
+def test_available_distances():
+    """available_distances has correct content"""
+    content = available_distances()
+    assert content.shape == (6, 2)
+    assert content["tn93", 1] == "dna, rna"
+
+
+def test_to_dict():
+    """distance matrix correctly produces a 1D dict"""
+    data = {("s1", "s2"): 0.25, ("s2", "s1"): 0.25}
+    dmat = DistanceMatrix(data)
+    got = dmat.to_dict()
+    assert got == data
+
+
+def test_matrix_dtype():
+    """tests DistanceMatrix correctly accepts the data with proper dtype"""
+    data = {
+        ("ABAYE2984", "Atu3667"): None,
+        ("ABAYE2984", "Avin_42730"): 0.638,
+        ("ABAYE2984", "BAA10469"): None,
+        ("Atu3667", "ABAYE2984"): None,
+        ("Atu3667", "Avin_42730"): 2.368,
+        ("Atu3667", "BAA10469"): None,
+        ("Avin_42730", "ABAYE2984"): 0.638,
+        ("Avin_42730", "Atu3667"): 2.368,
+        ("Avin_42730", "BAA10469"): 1.85,
+        ("BAA10469", "ABAYE2984"): None,
+        ("BAA10469", "Atu3667"): None,
+        ("BAA10469", "Avin_42730"): 1.85,
+    }
+    names = set()
+    for p in data:
+        names.update(p)
+
+    # tests when data has None values and DistanceMatrix using dtype('float')
+    darr = DistanceMatrix(data)
+    assert darr.shape == (4, 4)
+    assert set(darr.names) == names
+    for (a, b), dist in data.items():
+        if dist is None:
+            assert numpy.isnan(darr[a, b])
+        else:
+            assert_allclose(dist, darr[a, b])
+
+    data = {
+        ("ABAYE2984", "Atu3667"): "None",
+        ("ABAYE2984", "Avin_42730"): 0.638,
+        ("ABAYE2984", "BAA10469"): None,
+        ("Atu3667", "ABAYE2984"): None,
+        ("Atu3667", "Avin_42730"): 2.368,
+        ("Atu3667", "BAA10469"): "None",
+        ("Avin_42730", "ABAYE2984"): 0.638,
+        ("Avin_42730", "Atu3667"): 2.368,
+        ("Avin_42730", "BAA10469"): 1.85,
+        ("BAA10469", "ABAYE2984"): None,
+        ("BAA10469", "Atu3667"): None,
+        ("BAA10469", "Avin_42730"): 1.85,
+    }
+
+    # tests when data has str values and DistanceMatrix using dtype('float')
+    with pytest.raises(ValueError):
+        _ = DistanceMatrix(data)
+
+
+def test_EstimateDistances(tmp_path, al):
+    """testing (well, exercising at least), EstimateDistances"""
+    d = EstimateDistances(al, JC69())
+    d.run(show_progress=False)
+    canned_result = {
+        ("b", "e"): 0.440840,
+        ("c", "e"): 0.440840,
+        ("a", "c"): 0.088337,
+        ("a", "b"): 0.188486,
+        ("a", "e"): 0.440840,
+        ("b", "c"): 0.0883373,
+    }
+    result = d.get_pairwise_distances().to_dict()
+    assert_dists_approx_equal(canned_result, result)
+
+    # excercise writing to file
+    d.write(tmp_path / "junk.txt")
+
+
+def test_take_dists():
+    """subsets the distance matrix"""
+    data = {
+        ("ABAYE2984", "Atu3667"): 0.25,
+        ("ABAYE2984", "Avin_42730"): 0.638,
+        ("ABAYE2984", "BAA10469"): None,
+        ("Atu3667", "ABAYE2984"): 0.25,
+        ("Atu3667", "Avin_42730"): 2.368,
+        ("Atu3667", "BAA10469"): 0.25,
+        ("Avin_42730", "ABAYE2984"): 0.638,
+        ("Avin_42730", "Atu3667"): 2.368,
+        ("Avin_42730", "BAA10469"): 1.85,
+        ("BAA10469", "ABAYE2984"): 0.25,
+        ("BAA10469", "Atu3667"): 0.25,
+        ("BAA10469", "Avin_42730"): 1.85,
+    }
+    darr = DistanceMatrix(data)
+    got1 = darr.take_dists(["ABAYE2984", "Atu3667", "Avin_42730"])
+    got2 = darr.take_dists("BAA10469", negate=True)
+    assert_allclose(got1.array.astype(float), got2.array.astype(float))
+
+
+def test_build_phylogeny():
+    """build a NJ tree"""
+    from cogent3 import make_tree
+
+    dists = {
+        ("DogFaced", "FlyingFox"): 0.05,
+        ("DogFaced", "FreeTaile"): 0.14,
+        ("DogFaced", "LittleBro"): 0.16,
+        ("DogFaced", "TombBat"): 0.15,
+        ("FlyingFox", "DogFaced"): 0.05,
+        ("FlyingFox", "FreeTaile"): 0.12,
+        ("FlyingFox", "LittleBro"): 0.13,
+        ("FlyingFox", "TombBat"): 0.14,
+        ("FreeTaile", "DogFaced"): 0.14,
+        ("FreeTaile", "FlyingFox"): 0.12,
+        ("FreeTaile", "LittleBro"): 0.09,
+        ("FreeTaile", "TombBat"): 0.1,
+        ("LittleBro", "DogFaced"): 0.16,
+        ("LittleBro", "FlyingFox"): 0.13,
+        ("LittleBro", "FreeTaile"): 0.09,
+        ("LittleBro", "TombBat"): 0.12,
+        ("TombBat", "DogFaced"): 0.15,
+        ("TombBat", "FlyingFox"): 0.14,
+        ("TombBat", "FreeTaile"): 0.1,
+        ("TombBat", "LittleBro"): 0.12,
+    }
+    dists = DistanceMatrix(dists)
+    got = dists.quick_tree(show_progress=False)
+    expect = make_tree(
+        treestring="((TombBat,(DogFaced,FlyingFox)),LittleBro,FreeTaile)"
+    )
+    assert expect.same_topology(got)
+
+
+def test_names():
+    """names property works"""
+    data = {
+        ("ABAYE2984", "Atu3667"): 0.25,
+        ("ABAYE2984", "Avin_42730"): 0.638,
+        ("ABAYE2984", "BAA10469"): None,
+        ("Atu3667", "ABAYE2984"): 0.25,
+        ("Atu3667", "Avin_42730"): 2.368,
+        ("Atu3667", "BAA10469"): 0.25,
+        ("Avin_42730", "ABAYE2984"): 0.638,
+        ("Avin_42730", "Atu3667"): 2.368,
+        ("Avin_42730", "BAA10469"): 1.85,
+        ("BAA10469", "ABAYE2984"): 0.25,
+        ("BAA10469", "Atu3667"): 0.25,
+        ("BAA10469", "Avin_42730"): 1.85,
+    }
+    names = set()
+    for p in data:
+        names.update(p)
+    darr = DistanceMatrix(data)
+    assert set(darr.names) == names
+    darr = darr.drop_invalid()
+    for n in ("ABAYE2984", "BAA10469"):
+        names.remove(n)
+    assert set(darr.names) == names
+
+
+def test_EstimateDistancesWithMotifProbs(al):
+    """EstimateDistances with supplied motif probs"""
+    motif_probs = {"A": 0.1, "C": 0.2, "G": 0.2, "T": 0.5}
+    d = EstimateDistances(al, HKY85(), motif_probs=motif_probs)
+    d.run(show_progress=False)
+    canned_result = {
+        ("a", "c"): 0.07537,
+        ("b", "c"): 0.07537,
+        ("a", "e"): 0.39921,
+        ("a", "b"): 0.15096,
+        ("b", "e"): 0.39921,
+        ("c", "e"): 0.37243,
+    }
+    result = d.get_pairwise_distances().to_dict()
+    assert_dists_approx_equal(canned_result, result, atol=1e-5)
+
+
+def test_EstimateDistances_fromThreeway(al):
+    """testing (well, exercising at least), EsimateDistances fromThreeway"""
+    d = EstimateDistances(al, JC69(), threeway=True)
+    d.run(show_progress=False)
+    canned_result = {
+        ("b", "e"): 0.495312,
+        ("c", "e"): 0.479380,
+        ("a", "c"): 0.089934,
+        ("a", "b"): 0.190021,
+        ("a", "e"): 0.495305,
+        ("b", "c"): 0.0899339,
+    }
+    result = d.get_pairwise_distances(summary_function="mean").to_dict()
+    assert_dists_approx_equal(canned_result, result, atol=1e-4)
+
+
+def test_EstimateDistances_fromUnaligned(collection):
+    """Excercising estimate distances from unaligned sequences"""
+    d = EstimateDistances(collection, JC69(), do_pair_align=True, rigorous_align=True)
+    d.run(show_progress=False)
+    canned_result = {
+        ("b", "e"): 0.440840,
+        ("c", "e"): 0.440840,
+        ("a", "c"): 0.088337,
+        ("a", "b"): 0.188486,
+        ("a", "e"): 0.440840,
+        ("b", "c"): 0.0883373,
+    }
+    result = d.get_pairwise_distances().to_dict()
+    assert_dists_approx_equal(canned_result, result)
+
+    d = EstimateDistances(collection, JC69(), do_pair_align=True, rigorous_align=False)
+    d.run(show_progress=False)
+    canned_result = {
+        ("b", "e"): 0.440840,
+        ("c", "e"): 0.440840,
+        ("a", "c"): 0.088337,
+        ("a", "b"): 0.188486,
+        ("a", "e"): 0.440840,
+        ("b", "c"): 0.0883373,
+    }
+    result = d.get_pairwise_distances().to_dict()
+    assert_dists_approx_equal(canned_result, result)
+
+
+def test_EstimateDistances_other_model_params(al):
+    """test getting other model params from EstimateDistances"""
+    d = EstimateDistances(al, HKY85(), est_params=["kappa"])
+    d.run(show_progress=False)
+    # this will be a Number object with Mean, Median etc ..
+    kappa = d.get_param_values("kappa")
+    assert_allclose(kappa.mean, 0.8939, atol=1e-4)
+    # this will be a dict with pairwise instances, it's called by the above
+    # method, so the correctness of it's values is already checked
+    _ = d.get_pairwise_param("kappa")
+
+
+def test_EstimateDistances_modify_lf(al):
+    """tests modifying the lf"""
+
+    def constrain_fit(lf):
+        lf.set_param_rule("kappa", is_constant=True)
+        lf.optimise(local=True, show_progress=False)
+        return lf
+
+    d = EstimateDistances(al, HKY85(), modify_lf=constrain_fit)
+    d.run(show_progress=False)
+    result = d.get_pairwise_distances().to_dict()
+    d = EstimateDistances(al, F81())
+    d.run(show_progress=False)
+    expect = d.get_pairwise_distances().to_dict()
+    assert_dists_approx_equal(expect, result)
+
+
+def test_get_raw_estimates(al):
+    """correctly return raw result object"""
+    d = EstimateDistances(al, HKY85(), est_params=["kappa"])
+    d.run(show_progress=False)
+    expect = {
+        ("a", "b"): {
+            "kappa": 1.0000226766004808e-06,
+            "length": 0.18232155856115662,
+        },
+        ("a", "c"): {
+            "kappa": 1.0010380037049357e-06,
+            "length": 0.087070406623635604,
+        },
+        ("a", "e"): {"kappa": 2.3965871843412687, "length": 0.4389176272584539},
+        ("b", "e"): {"kappa": 2.3965871854366592, "length": 0.43891762729173389},
+        ("b", "c"): {
+            "kappa": 1.0010380037049357e-06,
+            "length": 0.087070406623635604,
+        },
+        ("c", "e"): {"kappa": 0.57046787478038707, "length": 0.43260232210282784},
+    }
+    got = d.get_all_param_values()
+    for pair in expect:
+        for param in expect[pair]:
+            assert_allclose(got[pair][param], expect[pair][param], atol=1e-6)
+
+
+def test_no_calc():
+    """returns None if no calculation done"""
+    al = load_aligned_seqs("data/brca1_5.paml")
+    d = EstimateDistances(al, submodel=HKY85())
+    assert d.get_pairwise_distances() is None
+
+
+def test_to_table():
+    """converts a distance matrix to a Table"""
+    data = {
+        ("A", "B"): 2,
+        ("A", "C"): 3,
+        ("B", "C"): 1,
+        ("B", "A"): 2,
+        ("C", "A"): 3,
+        ("C", "B"): 1,
+    }
+    darr = DistanceMatrix(data)
+    table = darr.to_table()
+    assert table.shape == (3, 4)
+    assert table.columns["names"].tolist() == list(darr.names)
+    assert table["A", "B"] == 2
+    assert table["A", "A"] == 0
+
+
+@pytest.mark.parametrize(
+    "aln", ["basic_alignment", "ambig_alignment", "diff_alignment"]
+)
+def test_jc69_dists(request, aln):
+    """full numba implementation matches original"""
+    aln = request.getfixturevalue(aln)
+    calc = JC69Pair(DNA, alignment=aln)
+    calc.run(show_progress=False)
+    expect = calc.dists.array
+    got = pdist_numba.jc69(aln)
+    assert_allclose(expect, got)
+
+
+@pytest.mark.parametrize("aln", ["basic_alignment", "diff_alignment"])
+def test_tn93_dists(request, aln):
+    """full numba implementation matches original"""
+    # note: I excluded the ambig_alignment case because the
+    # old and new approaches differ in how they compute the
+    # nucleotide frequencies. The old approach would effectively
+    # exclude data in aligned columns that contained a
+    # non-canonical state, the approach computes the nucleotide
+    # frequencies from the whole alignment and so does not exclude
+    # valid states just because another sequence has an ambiguity
+    aln = request.getfixturevalue(aln)
+    # convert to new type alignment
+    aln = make_aligned_seqs(data=aln.to_dict(), moltype="dna", new_type=True)
+    calc = TN93Pair(aln.moltype, alignment=aln)
+    calc.run(show_progress=False)
+    expect = calc.dists.array
+    got = pdist_numba.tn93(aln, parallel=False)
+    assert_allclose(got.array, expect)
+
+
+@pytest.mark.parametrize("aln", ["basic_alignment", "diff_alignment"])
+def test_paralinear_dists(request, aln):
+    aln = request.getfixturevalue(aln)
+    # convert to new type alignment
+    aln = make_aligned_seqs(data=aln.to_dict(), moltype="dna", new_type=True)
+    calc = ParalinearPair(aln.moltype, alignment=aln)
+    calc.run(show_progress=False)
+    expect = calc.dists.array
+    got = pdist_numba.paralinear(aln, parallel=False)
+    assert_allclose(got.array, expect, atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    "aln", ["basic_alignment", "ambig_alignment", "diff_alignment"]
+)
+def test_hamming_dists(request, aln):
+    aln = request.getfixturevalue(aln)
+    # convert to new type alignment
+    aln = make_aligned_seqs(data=aln.to_dict(), moltype="dna", new_type=True)
+    calc = HammingPair(aln.moltype, alignment=aln)
+    calc.run(show_progress=False)
+    expect = calc.dists.array
+    got = pdist_numba.hamming(aln, parallel=False)
+    assert_allclose(got.array, expect)
+
+
+@pytest.mark.parametrize(
+    "aln", ["basic_alignment", "ambig_alignment", "diff_alignment"]
+)
+def test_prop_dists(request, aln):
+    aln = request.getfixturevalue(aln)
+    # convert to new type alignment
+    aln = make_aligned_seqs(data=aln.to_dict(), moltype="dna", new_type=True)
+    calc = ProportionIdenticalPair(aln.moltype, alignment=aln)
+    calc.run(show_progress=False)
+    expect = calc.dists.array
+    got = pdist_numba.pdist(aln, parallel=False)
+    assert_allclose(got.array, expect)
+
+
+@pytest.mark.parametrize("calc", ["jc69", "tn93", "paralinear", "hamming", "pdist"])
+def test_numba_get_dist(calc, basic_alignment) -> None:
+    aln = make_aligned_seqs(
+        data=basic_alignment.to_dict(), moltype="dna", new_type=True
+    )
+    calc = pdist_numba.get_distance_calculator(calc)
+    got = calc(aln)
+    assert isinstance(got, DistanceMatrix)
+
+
+@pytest.mark.parametrize("calc", ["tn93", "paralinear"])
+def test_invalid_moltype_fast_distances(calc, basic_alignment):
+    from cogent3.core import new_moltype
+
+    aln = make_aligned_seqs(
+        data=basic_alignment.to_dict(), moltype="protein", new_type=True
+    )
+    calc = pdist_numba.get_distance_calculator(calc)
+    with pytest.raises(new_moltype.MolTypeError):
+        calc(aln, invalid_raises=True)
+
+
+def test_new_paralinear_for_determinant_lte_zero():
+    """returns distance of None if the determinant is <= 0"""
+    data = dict(
+        seq1="AGGGGGGGGGGCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGCGGTTTTTTTTTTTTTTTTTT",
+        seq2="TAAAAAAAAAAGGGGGGGGGGGGGGGGGGTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCC",
+    )
+    aln = make_aligned_seqs(data=data, moltype="dna", new_type=True)
+
+    dists = pdist_numba.paralinear(aln, parallel=False)
+    assert numpy.isnan(dists.array).any()
+
+    # raises ArithmeticError if told to
+    with pytest.raises(ArithmeticError):
+        dists = pdist_numba.paralinear(aln, parallel=False, invalid_raises=True)
+
+    # also raised if alignment method used
+    with pytest.raises(ArithmeticError):
+        aln.distance_matrix(calc="paralinear", drop_invalid=False)
+
+
+def test_unknown_dist():
+    with pytest.raises(ValueError):
+        pdist_numba.get_distance_calculator("blah")
+
+
+@pytest.mark.parametrize("calc", ["jc69", "tn93", "paralinear", "hamming", "pdist"])
+def test_compare_parallel_serial(DATA_DIR, calc):
+    aln = load_aligned_seqs(DATA_DIR / "brca1.fasta", moltype="dna", new_type=True)
+    serial = aln.distance_matrix(calc=calc, parallel=False)
+    parallel = aln.distance_matrix(calc=calc, parallel=True)
+    assert_allclose(serial.array, parallel.array)