Repository: cogent3/cogent3
Pull request #1808 (2024-04-08 04:26:18)

MAINT: merged changes on develop branch into map refactor


Patch:
diff --git a/.github/ISSUE_TEMPLATE/deprecation.md b/.github/ISSUE_TEMPLATE/deprecation.md
index dd774623e..5fa83d045 100644
--- a/.github/ISSUE_TEMPLATE/deprecation.md
+++ b/.github/ISSUE_TEMPLATE/deprecation.md
@@ -21,9 +21,10 @@ What should it be replaced with and why?
 
 **How much notice**
 
-The number of months (*x*) long before the capability will be removed? This should be represented in the user warning as a specific version. e.g. <current year>.<current month + *x*> where *x* is from the PR date.  
+**How much notice**
 
-**Guidelines for deprecation**
+How much time should the community be given to response to this notice? Specify as the number of months (*x*) before the capability will be removed. Leave blank if not known.
 
-As this is a function being deprecated, see [dev guidelines for deprecation](https://github.com/cogent3/c3dev/wiki/Deprecating-Code) of code. Indicate the limiting version (for when the old function will be removed) as 3 months from when you make the change.
+**Guidelines for deprecation**
 
+As this is a function that will be removed, see [dev guidelines for deprecation](https://github.com/cogent3/c3dev/wiki/Deprecating-Code) of code. Indicate the limiting version (for when the old function will be removed) as <current year>.<PR month + *x*> where *x* is defined above.
diff --git a/.github/ISSUE_TEMPLATE/discontinue.md b/.github/ISSUE_TEMPLATE/discontinue.md
index 28740a010..8f82f7e26 100644
--- a/.github/ISSUE_TEMPLATE/discontinue.md
+++ b/.github/ISSUE_TEMPLATE/discontinue.md
@@ -21,9 +21,8 @@ Why should it be discontinued?
 
 **How much notice**
 
-The number of months (*x*) long before the capability will be removed? This should be represented in the user warning as a specific version. e.g. <current year>.<current month + *x*> where *x* is from the PR date.  
+How much time should the community be given to response to this notice? Specify as the number of months (*x*) before the capability will be removed. Leave blank if not known.
 
 **Guidelines for deprecation**
 
-As this is a function that will be removed, see [dev guidelines for deprecation](https://github.com/cogent3/c3dev/wiki/Deprecating-Code) of code. Indicate the limiting version (for when the old function will be removed) as 3 months from when you make the change.
-
+As this is a function that will be removed, see [dev guidelines for deprecation](https://github.com/cogent3/c3dev/wiki/Deprecating-Code) of code. Indicate the limiting version (for when the old function will be removed) as <current year>.<PR month + *x*> where *x* is defined above.  
diff --git a/.github/workflows/linters.yml b/.github/workflows/linters.yml
index 81ac7f97a..432824952 100644
--- a/.github/workflows/linters.yml
+++ b/.github/workflows/linters.yml
@@ -12,7 +12,7 @@ on:
 jobs:
   linters:
     runs-on: ubuntu-latest
-    if: github.repository == 'cogent/cogent3'
+    if: github.repository == 'cogent3/cogent3'
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python
@@ -33,7 +33,7 @@ jobs:
       - name: Format code using black
         run: black .
       - name: Sort import statements using isort
-        run: isort
+        run: isort .
       - name: Commit changes
         uses: EndBug/add-and-commit@v9
         with:
diff --git a/.gitignore b/.gitignore
index f557c8489..3bd13a6e5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,6 +19,9 @@ lib
 lib64
 doc/_build
 
+# docker stuff
+.devcontainer/*
+
 # Installer logs
 pip-log.txt
 
diff --git a/.hgignore b/.hgignore
index b5a6403a3..76a64dda6 100644
--- a/.hgignore
+++ b/.hgignore
@@ -37,3 +37,4 @@ dist/*
 working/*
 lcov*.info
 .ruff_cache/*
+.devcontainer/*
\ No newline at end of file
diff --git a/changelog.d/20240325_161732_Gavin.Huttley.md b/changelog.d/20240325_161732_Gavin.Huttley.md
new file mode 100644
index 000000000..65d9410b9
--- /dev/null
+++ b/changelog.d/20240325_161732_Gavin.Huttley.md
@@ -0,0 +1,43 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+-->
+
+### Contributors
+
+- First-time contributions from Yutong-Shao, Zongjing-Han, yantonglu,
+  ShunyiYang, KyleSu12, Ruizhe-wang, Firimp 🎉
+- Adopted suggested documentation example proposed by Lizhaozhe123,
+  with some amendments, for the trim_stop_codons app.
+
+<!--
+### ENH
+
+- A bullet item for the ENH category.
+
+-->
+<!--
+### BUG
+
+- A bullet item for the BUG category.
+
+-->
+<!--
+### DOC
+
+- A bullet item for the DOC category.
+
+-->
+<!--
+### Deprecations
+
+- A bullet item for the Deprecations category.
+
+-->
+<!--
+### Discontinued
+
+- A bullet item for the Discontinued category.
+
+-->
diff --git a/changelog.d/20240328_134342_Gavin.Huttley.md b/changelog.d/20240328_134342_Gavin.Huttley.md
new file mode 100644
index 000000000..9462e2619
--- /dev/null
+++ b/changelog.d/20240328_134342_Gavin.Huttley.md
@@ -0,0 +1,47 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+-->
+
+
+### Contributors
+
+- rmcar17 identified and fixed a bug 🚀!
+- rmcar17 addressed deprecation of using `end` as a column name
+
+
+<!--
+### ENH
+
+- A bullet item for the ENH category.
+
+-->
+
+### BUG
+
+- The saving of numpy arrays in annotation db's was not cross-platform
+  compatible due to differences in default types between OS's. Fixed by
+  using numpy.save and numpy.load. A warning is raised if an old style
+  format is detected and a function provided for updating the format.
+
+
+<!--
+### DOC
+
+- A bullet item for the DOC category.
+
+-->
+
+### Deprecations
+
+- `end` is a SQL keyword and so should be avoided. It is now replaced
+  by `stop`. Annotation db's with the old column name are dynamically
+  updated.
+
+<!--
+### Discontinued
+
+- A bullet item for the Discontinued category.
+
+-->
diff --git a/changelog.d/20240329_152534_Gavin.Huttley.md b/changelog.d/20240329_152534_Gavin.Huttley.md
new file mode 100644
index 000000000..779132349
--- /dev/null
+++ b/changelog.d/20240329_152534_Gavin.Huttley.md
@@ -0,0 +1,48 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+-->
+
+<!--
+### Contributors
+
+- A bullet item for the Contributors category.
+
+-->
+
+### ENH
+
+- methods on likelihood function objects that work with continuous-time
+  Markov models no longer fail if the model also has discrete-time edges.
+  These include lf.get_lengths_as_ens(), lf.get_annotated_tree(),
+  lf.get_paralinear_metric().
+- Added new lf.get_ens_tree(). This returns trees with the expected number
+  of substitutions as the branch length. On discrete-time edges the branch
+  length is set to None. Think of this tree as the true "evolutionary tree".
+  Thanks to Von Bing Yap for suggesting this!
+
+<!--
+### BUG
+
+- A bullet item for the BUG category.
+
+-->
+<!--
+### DOC
+
+- A bullet item for the DOC category.
+
+-->
+<!--
+### Deprecations
+
+- A bullet item for the Deprecations category.
+
+-->
+<!--
+### Discontinued
+
+- A bullet item for the Discontinued category.
+
+-->
diff --git a/requirements.txt b/requirements.txt
index 18ab6e666..03b4568a2 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,7 +7,7 @@ numba>0.53; python_version>='3.9'
 numba>0.54; python_version>='3.10'
 numba>=0.57.0; python_version=='3.11'
 scipy==1.12.0
-nox==2023.4.22
-plotly==5.19.0
+nox==2024.3.2
+plotly==5.20.0
 typing_extensions>=3.7
 .[dev]
diff --git a/src/cogent3/app/composable.py b/src/cogent3/app/composable.py
index 0496508ae..9bab28057 100644
--- a/src/cogent3/app/composable.py
+++ b/src/cogent3/app/composable.py
@@ -788,7 +788,7 @@ def _apply_to(
         inputs, parallel=parallel, par_kw=par_kw, show_progress=show_progress
     ):
         member = self.main(data=result.obj, identifier=id_from_source(result.source))
-        md5 = member.md5
+        md5 = getattr(member, "md5", None)
         logger.log_message(str(member), label="output")
         if md5:
             logger.log_message(md5, label="output md5sum")
diff --git a/src/cogent3/app/result.py b/src/cogent3/app/result.py
index a2d23cbc8..aacdc75e1 100644
--- a/src/cogent3/app/result.py
+++ b/src/cogent3/app/result.py
@@ -341,6 +341,7 @@ class model_result(generic_result):
     @property
     def tree(self):
         """an annotated tree with 'ENS' set as the branch length
+
         Note
         ----
         In the case of a discrete time process, length is 'paralinear'"""
diff --git a/src/cogent3/app/sample.py b/src/cogent3/app/sample.py
index e47ad759d..1dccf4b2e 100644
--- a/src/cogent3/app/sample.py
+++ b/src/cogent3/app/sample.py
@@ -327,24 +327,66 @@ class take_codon_positions:
 
     def __init__(
         self,
-        *positions,
-        fourfold_degenerate=False,
-        gc="Standard",
-        moltype="dna",
+        *positions: int,
+        fourfold_degenerate: bool = False,
+        gc: Union[str, int] = "Standard",
+        moltype: str = "dna",
     ):
         """
         Parameters
         ----------
         positions
-            either an integer (1, 2, 3), or a tuple of position numbers,
-            e.g. 3 is third position, (1,2) is first and second codon position
-        fourfold_degenerate : bool
+            either a single integer from (1, 2, 3), or additional keyword
+            arguments of position numbers, e.g. 3 is third position, (1,2)
+            is first and second codon position
+        fourfold_degenerate
             if True, returns third positions from four-fold degenerate codons.
             Overrides positions.
         gc
-            identifier for a genetic code or a genetic code instance
-        moltype : str
+            identifier for a genetic code or a genetic code instance.
+            see https://cogent3.org/doc/cookbook/what_codes.html
+        moltype
             molecular type, must be either DNA or RNA
+
+        Examples
+        --------
+
+        Create a sample alignment and an app that extracts the 3rd codon
+        position from an alignment.
+
+        >>> from cogent3 import make_aligned_seqs, get_app
+        >>> aln = make_aligned_seqs({"s1": "ACGACGACG", "s2": "GATGATGAT"})
+        >>> take_pos3 = get_app("take_codon_positions", 3, moltype="dna")
+        >>> result = take_pos3(aln)
+        >>> print(result.to_pretty())
+        s1    GGG
+        s2    TTT
+
+        Create an app that extracts the 1st and 2nd codon positions from an
+        alignment.
+
+        >>> take_pos12 = get_app("take_codon_positions", 1, 2, moltype="dna")
+        >>> result = take_pos12(aln)
+        >>> print(result.to_pretty())
+        s1    ACACAC
+        s2    GAGAGA
+
+        Create a sample alignment and an app that returns the 3rd codon
+        positions from four-fold degenerate codons.
+
+        >>> aln_ff = make_aligned_seqs({"s1": "GCAAGCGTTTAT", "s2": "GCTTTTGTCAAT"})
+        >>> take_fourfold = get_app("take_codon_positions", fourfold_degenerate=True, moltype="dna")
+        >>> result = take_fourfold(aln_ff)
+        >>> print(result.to_pretty())
+        s1    AT
+        s2    TC
+
+        A NotCompleted object (see https://cogent3.org/doc/app/not-completed.html)
+        is returned if all sites are excluded.
+
+        >>> result = take_fourfold(aln)
+        >>> result.message
+        'Traceback ...
         """
         assert moltype is not None
         moltype = get_moltype(moltype)
@@ -754,6 +796,61 @@ class omit_duplicated:
             set random number seed. Only applied of choose=='random'
         moltype
             molecular type, can be string or instance
+
+        Examples
+        --------
+        Removes redundant sequences from a sequence collection (aligned or
+        unaligned).
+
+        Create sample data with duplicated sequences.
+
+        >>> from cogent3 import get_app, make_aligned_seqs
+        >>> data = {
+        ... "a": "ACGT",
+        ... "b": "ACG-",  # identical to 'a' except has a gap
+        ... "c": "ACGG",  # duplicate
+        ... "d": "ACGG",  # duplicate
+        ... "e": "AGTC"   # unique
+        ... }
+
+        Create an app that selects a representative of omits duplicate sequences.
+        Setting ``choose="longest"`` selects the duplicated sequence with the least
+        number of gaps and ambiguous characters. In this case, only one of 'c' and
+        'd' will be retained.
+
+        >>> seqs = make_aligned_seqs(data, moltype="DNA")
+        >>> app = get_app("omit_duplicated", moltype="dna", choose="longest")
+        >>> result = app(seqs)
+        >>> print(result.to_pretty())
+        a    ACGT
+        b    ...-
+        d    ...G
+        e    .GTC
+
+        Setting ``choose=None`` means only unique sequences are retained.
+
+        >>> app = get_app("omit_duplicated", moltype="dna", choose=None)
+        >>> result = app(seqs)
+        >>> print(result.to_pretty())
+        a    ACGT
+        b    ...-
+        e    .GTC
+
+        Use the ``mask_degen`` argument to specify how to treat matches between
+        sequences with degenerate characters. We create sample data first that
+        has a DNA ambiguity code.
+
+        >>> data = make_aligned_seqs({
+        ... "s1": "ATCG",
+        ... "s2": "ATYG",  # matches s1 with ambiguity
+        ... "s3": "GGTA",
+        ... },
+        ... moltype="DNA")
+        >>> app_dna = get_app("omit_duplicated", mask_degen=True, choose="longest", moltype="DNA")
+        >>> result = app_dna(data)
+        >>> print(result.to_pretty())
+        s1    ATCG
+        s3    GGTA
         """
         assert not choose or choose in "longestrandom"
         if moltype:
@@ -820,18 +917,46 @@ class omit_duplicated:
 class trim_stop_codons:
     """Removes terminal stop codons."""
 
-    def __init__(self, gc=1):
+    def __init__(self, gc: Union[str, int] = 1):
         """
         Parameters
         ----------
         gc
             identifier for a genetic code or a genetic code instance, defaults
-            to standard genetic code
+            to standard genetic code. See https://cogent3.org/doc/cookbook/what_codes.html
 
         Returns
         -------
         A new sequence collection, or False if not all the named sequences are
         in the collection.
+
+        Examples
+        --------
+
+        Removes trailing (terminal) stop codons from sequences in a sequence
+        collection or alignment.
+
+        Create a sample unaligned sequence data and an app to remove trailing
+        stop codons.
+
+        >>> from cogent3 import make_unaligned_seqs, make_aligned_seqs, get_app
+        >>> ualn = make_unaligned_seqs(data={"s1": "AAATTTCCC", "s2": "AAATTTTAA"}, moltype="dna")
+        >>> app = get_app("trim_stop_codons")
+        >>> result = app(ualn)
+        >>> print(result.to_dict())
+        {'s1': 'AAATTTCCC', 's2': 'AAATTT'}
+
+        Remove trailing stop codons from aligned sequences, using the genetic
+        code for "Vertebrate Mitochondrial".
+
+        For a list of all genetic codes, see
+        https://cogent3.org/doc/cookbook/what_codes.html.
+
+        >>> aln = make_aligned_seqs(data={"s1": "AAATTTCCC", "s2": "AAATTTTAA"}, moltype="dna")
+        >>> app = get_app("trim_stop_codons", gc=2)
+        >>> result = app(aln)
+        >>> print(result.to_dict())
+        {'s1': 'AAATTTCCC', 's2': 'AAATTT---'}
         """
         self._gc = gc
 
diff --git a/src/cogent3/app/typing.py b/src/cogent3/app/typing.py
index 660fdd076..ca2c7bc2b 100644
--- a/src/cogent3/app/typing.py
+++ b/src/cogent3/app/typing.py
@@ -42,7 +42,6 @@ ResultType = Union[
     TabularResultType,
 ]
 
-# todo when move to python 3.8 define protocols for IdentifierType and SerialisableType
 IdentifierType = TypeVar("IdentifierType")
 
 
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 76f0c293f..3ba2a0f6a 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -18,6 +18,7 @@
     passed in a stream of two-item label, sequence pairs. However, this can
     cause confusion when testing.
 """
+
 from __future__ import annotations
 
 import functools
@@ -986,15 +987,7 @@ class _SequenceCollectionBase:
         return json.dumps(self.to_rich_dict())
 
     def to_fasta(self):
-        """Return alignment in Fasta format
-
-        Parameters
-        ----------
-        make_seqlabel
-            callback function that takes the seq object and
-            returns a label str
-
-        """
+        """Return alignment in Fasta format"""
         return alignment_to_fasta(self.to_dict())
 
     def to_nexus(self, seq_type, wrap=50):
@@ -3936,11 +3929,11 @@ class ArrayAlignment(AlignmentI, _SequenceCollectionBase):
         self._seqs = curr_seqs
         self.seq_len = curr_seqs.shape[1] if len(curr_seqs) else 0
 
-    def _get_positions(self):
+    @property
+    def positions(self):
         """Override superclass positions to return positions as symbols."""
-        return list(map(self.alphabet.from_indices, self.array_positions))
-
-    positions = property(_get_positions)
+        from_indices = self.alphabet.from_indices
+        return [list(from_indices(pos)) for pos in self.array_positions]
 
     @property
     def named_seqs(self):
@@ -5295,7 +5288,7 @@ class Alignment(AlignmentI, SequenceCollection):
             seqname = seqid_to_seqname[seqid]
             seq = self.named_seqs[seqname]
             # we use parent seqid, stored on SeqView
-            parent_id, start, end, _ = seq.data.parent_coordinates()
+            parent_id, start, stop, _ = seq.data.parent_coordinates()
             offset = seq.data.annotation_offset
 
             for feature in self.annotation_db.get_features_matching(
@@ -5305,7 +5298,7 @@ class Alignment(AlignmentI, SequenceCollection):
                 on_alignment=False,
                 allow_partial=allow_partial,
                 start=start,
-                end=end,
+                stop=stop,
             ):
                 if offset:
                     feature["spans"] = (array(feature["spans"]) - offset).tolist()
@@ -5534,6 +5527,15 @@ def _coerce_to_array_aligned_seqs(data, names, label_to_name=str, moltype=None)
         seq = _construct_array_aligned_seq(data[name], moltype=moltype)
         seqs.append(seq)
 
+    # Ensure all sequences are of equal length
+    seq_lengths = {len(seq) for seq in seqs}
+    if len(seq_lengths) > 1:
+        raise ValueError(
+            f"Input sequences are not all the same length: {seq_lengths}. "
+            "Please ensure all sequences are properly aligned or "
+            "correct the input file format."
+        )
+
     return numpy.array(seqs), names
 
 
diff --git a/src/cogent3/core/alphabet.py b/src/cogent3/core/alphabet.py
index 89721c313..452520440 100644
--- a/src/cogent3/core/alphabet.py
+++ b/src/cogent3/core/alphabet.py
@@ -16,11 +16,13 @@ and/or Enumerations on the fly, however.
 import json
 import typing
 
+from functools import singledispatchmethod
 from itertools import product
 
 from numpy import (
     arange,
     array,
+    ndarray,
     newaxis,
     sum,
     take,
@@ -624,7 +626,7 @@ class CharAlphabet(Alphabet):
         data = data or []
         super(CharAlphabet, self).__init__(data, gap, moltype=moltype)
         self._indices_to_chars, self._chars_to_indices = _make_translation_tables(data)
-        self._char_nums_to_indices = array(range(256), uint8)
+        self._char_nums_to_indices = array(range(256), self.array_type)
         for c, i in self._chars_to_indices.items():
             self._char_nums_to_indices[c] = i
 
@@ -672,6 +674,76 @@ class CharAlphabet(Alphabet):
                 [i.tobytes().decode("utf-8") for i in self.to_chars(data)]
             )
 
+    @singledispatchmethod
+    def from_indices(self, data: ndarray) -> str:
+        """returns string from numpy array of integers"""
+        data = data.astype(dtype=self.array_type)
+        encoding = data.dtype.name.replace("uint", "utf")
+        return str.translate(data.tobytes().decode(encoding), self._indices_to_chars)
+
+    @from_indices.register
+    def _(self, data: list) -> str:
+        """returns string from numpy array of integers"""
+        return self.from_indices(array(data, dtype=self.array_type))
+
+    @from_indices.register
+    def _(self, data: tuple) -> str:
+        """returns string from numpy array of integers"""
+        return self.from_indices(array(data, dtype=self.array_type))
+
+    @singledispatchmethod
+    def to_indices(self, data) -> ndarray:
+        """Returns sequence of indices from sequence of elements.
+
+        Raises KeyError if some of the elements were not found.
+
+        Expects data to be a sequence (e.g. list of tuple) of items that
+        are in the Enumeration. Returns a list containing the index of each
+        element in the input, in order.
+
+        e.g. for the RNA alphabet ('U','C','A','G'), the sequence 'CCAU'
+        would produce the result [1,1,2,0], returning the index of each
+        element in the input.
+        """
+        return array([self._obj_to_index[e] for e in data], dtype=self.array_type)
+
+    @to_indices.register
+    def _(self, data: str) -> ndarray:
+        """Returns sequence of indices from sequence of elements.
+
+        Raises KeyError if some of the elements were not found.
+
+        Expects data to be a sequence (e.g. list of tuple) of items that
+        are in the Enumeration. Returns a list containing the index of each
+        element in the input, in order.
+
+        e.g. for the RNA alphabet ('U','C','A','G'), the sequence 'CCAU'
+        would produce the result [1,1,2,0], returning the index of each
+        element in the input.
+        """
+        return array(
+            memoryview(
+                bytearray(data.translate(self._chars_to_indices).encode("utf8"))
+            ),
+            dtype=self.array_type,
+        )
+
+    @to_indices.register
+    def _(self, data: bytes) -> ndarray:
+        """Returns sequence of indices from sequence of elements.
+
+        Raises KeyError if some of the elements were not found.
+
+        Expects data to be a sequence (e.g. list of tuple) of items that
+        are in the Enumeration. Returns a list containing the index of each
+        element in the input, in order.
+
+        e.g. for the RNA alphabet ('U','C','A','G'), the sequence 'CCAU'
+        would produce the result [1,1,2,0], returning the index of each
+        element in the input.
+        """
+        return self.to_indices(data.decode("utf8"))
+
 
 T = typing.Union[Alphabet, CharAlphabet]
 
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index 0c70f8263..fa2b8f8c2 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -171,7 +171,7 @@ class Feature:
         """generator returns sub-features of self optionally matching biotype"""
         offset = getattr(self.parent, "annotation_offset", 0)
         start = self.map.start + offset
-        end = self.map.end + offset
+        stop = self.map.end + offset
 
         make_feature = self.parent.make_feature
         db = self.parent.annotation_db
@@ -179,7 +179,7 @@ class Feature:
             biotype=biotype,
             name=self.name,
             start=start,
-            end=end,
+            stop=stop,
             **kwargs,
         ):
             record["spans"] = array(record["spans"]) - offset
@@ -195,14 +195,14 @@ class Feature:
         """generator returns parent features of self optionally matching biotype"""
         offset = getattr(self.parent, "annotation_offset", 0)
         start = self.map.start + offset
-        end = self.map.end + offset
+        stop = self.map.end + offset
 
         make_feature = self.parent.make_feature
         db = self.parent.annotation_db
         for record in db.get_feature_parent(
             name=self.name,
             start=start,
-            end=end,
+            stop=stop,
             **kwargs,
         ):
             record["spans"] = array(record["spans"]) - offset
diff --git a/src/cogent3/core/annotation_db.py b/src/cogent3/core/annotation_db.py
index 9e05b2e7f..60bd064ed 100644
--- a/src/cogent3/core/annotation_db.py
+++ b/src/cogent3/core/annotation_db.py
@@ -3,6 +3,7 @@ from __future__ import annotations
 import collections
 import copy
 import inspect
+import io
 import json
 import os
 import pathlib
@@ -10,10 +11,13 @@ import re
 import sqlite3
 import sys
 import typing
+import warnings
 
 import numpy
 import typing_extensions
 
+import cogent3.util.warning as c3warn
+
 from cogent3._version import __version__
 from cogent3.util.deserialise import register_deserialiser
 from cogent3.util.misc import extend_docstring_from, get_object_provenance
@@ -36,14 +40,35 @@ _is_ge_3_11 = (sys.version_info.major, sys.version_info.minor) >= (3, 11)
 
 # Define custom types for storage in sqlite
 # https://stackoverflow.com/questions/18621513/python-insert-numpy-array-into-sqlite3-database
-def array_to_sqlite(data):
-    return data.tobytes()
-
+def array_to_sqlite(data: numpy.ndarray) -> bytes:
+    with io.BytesIO() as out:
+        numpy.save(out, data)
+        out.seek(0)
+        output = out.read()
+    return output
+
+
+def sqlite_to_array(data: bytes) -> numpy.ndarray:
+    with io.BytesIO(data) as out:
+        out.seek(0)
+        try:
+            result = numpy.load(out)
+        except ValueError:
+            # array is not stored in the numpy.save format
+            # attempt to read from the old format where the
+            # array was saved using numpy.ndarray.tobytes
+            result = numpy.frombuffer(data, dtype=int)
+            dim = result.shape[0] // 2
+            result = result.reshape((dim, 2))
+
+            warnings.warn(
+                "Old OS dependent database file format detected. "
+                "Update the file format using cogent3.core.annotation_db.update_file_format() "
+                "For reason see https://github.com/cogent3/cogent3/issues/1776.",
+                UserWarning,
+            )
 
-def sqlite_to_array(data):
-    result = numpy.frombuffer(data, dtype=int)
-    dim = result.shape[0] // 2
-    return result.reshape((dim, 2))
+    return result
 
 
 def dict_to_sqlite_as_json(data: dict) -> str:
@@ -85,6 +110,11 @@ class SerialisableType(typing.Protocol):  # pragma: no cover
 @typing.runtime_checkable
 class SupportsQueryFeatures(typing.Protocol):  # pragma: no cover
     # should be defined centrally
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_features_matching(
         self,
         *,
@@ -92,20 +122,40 @@ class SupportsQueryFeatures(typing.Protocol):  # pragma: no cover
         biotype: OptionalStr = None,
         name: OptionalStr = None,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
         strand: OptionalStr = None,
         attributes: OptionalStr = None,
         on_alignment: OptionalBool = None,
     ) -> typing.Iterator[FeatureDataType]:
         ...
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_feature_children(
-        self, *, name: str, start: OptionalInt = None, end: OptionalInt = None, **kwargs
+        self,
+        *,
+        name: str,
+        start: OptionalInt = None,
+        stop: OptionalInt = None,
+        **kwargs,
     ) -> typing.List[FeatureDataType]:
         ...
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_feature_parent(
-        self, *, name: str, start: OptionalInt = None, end: OptionalInt = None, **kwargs
+        self,
+        *,
+        name: str,
+        start: OptionalInt = None,
+        stop: OptionalInt = None,
+        **kwargs,
     ) -> typing.List[FeatureDataType]:
         ...
 
@@ -121,6 +171,11 @@ class SupportsQueryFeatures(typing.Protocol):  # pragma: no cover
     ) -> int:
         ...
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def subset(
         self,
         *,
@@ -128,7 +183,7 @@ class SupportsQueryFeatures(typing.Protocol):  # pragma: no cover
         biotype: OptionalStr = None,
         name: OptionalStr = None,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
         strand: OptionalStr = None,
         attributes: OptionalStr = None,
     ):
@@ -247,8 +302,8 @@ def _matching_conditions(
     conditions : dict
         column name and values to be matched
     allow_partial : bool, optional
-        if False, only records within start, end are included. If True,
-        all records that overlap the segment defined by start, end are included.
+        if False, only records within start, stop are included. If True,
+        all records that overlap the segment defined by start, stop are included.
 
     Returns
     -------
@@ -256,7 +311,7 @@ def _matching_conditions(
         the SQL statement and the tuple of values
     """
     start = conditions.pop("start", None)
-    end = conditions.pop("end", None)
+    stop = conditions.pop("stop", None)
 
     sql = []
     vals = ()
@@ -278,27 +333,27 @@ def _matching_conditions(
         sql.append(" AND ".join(conds))
         vals = tuple(vals)
 
-    if isinstance(start, int) and isinstance(end, int):
+    if isinstance(start, int) and isinstance(stop, int):
         if allow_partial:
             # allow matches that overlap the segment
             cond = [
-                f"(start >= {start} AND end <= {end})",  # lies within the segment
-                f"(start <= {start} AND end > {start})",  # straddles beginning of segment
-                f"(start < {end} AND end >= {end})",  # straddles end of segment
-                f"(start <= {start} AND end >= {end})",  # includes segment
+                f"(start >= {start} AND stop <= {stop})",  # lies within the segment
+                f"(start <= {start} AND stop > {start})",  # straddles beginning of segment
+                f"(start < {stop} AND stop >= {stop})",  # straddles stop of segment
+                f"(start <= {start} AND stop >= {stop})",  # includes segment
             ]
             cond = " OR ".join(cond)
         else:
             # only matches within bounds
-            cond = f"start >= {start} AND end <= {end}"
+            cond = f"start >= {start} AND stop <= {stop}"
         sql.append(f"({cond})")
     elif isinstance(start, int):
-        # if query has no end, then any feature containing start
-        cond = f"(start <= {start} AND {start} < end)"
+        # if query has no stop, then any feature containing start
+        cond = f"(start <= {start} AND {start} < stop)"
         sql.append(f"({cond})")
-    elif isinstance(end, int):
-        # if query has no start, then any feature containing end
-        cond = f"(start <= {end} AND {end} < end)"
+    elif isinstance(stop, int):
+        # if query has no start, then any feature containing stop
+        cond = f"(start <= {stop} AND {stop} < stop)"
         sql.append(f"({cond})")
 
     sql = f"{' AND '.join(sql)}"
@@ -309,7 +364,7 @@ def _del_records_sql(
     table_name: str,
     conditions: dict,
     start: OptionalInt = None,
-    end: OptionalInt = None,
+    stop: OptionalInt = None,
     allow_partial=True,
 ) -> ReturnType:
     """creates the SQL and values for identifying records to be deleted
@@ -320,12 +375,12 @@ def _del_records_sql(
         table to have records deleted from
     conditions : dict
         column name and values to be matched
-    start, end : OptionalInt
-        select records whose (start, end) values lie between start and end,
+    start, stop : OptionalInt
+        select records whose (start, stop) values lie between start and stop,
         or overlap them if (allow_partial is True)
     allow_partial : bool, optional
-        if False, only records within start, end are included. If True,
-        all records that overlap the segment defined by start, end are included.
+        if False, only records within start, stop are included. If True,
+        all records that overlap the segment defined by start, stop are included.
 
     Returns
     -------
@@ -333,7 +388,7 @@ def _del_records_sql(
         the SQL statement and the tuple of values
     """
     where, vals = _matching_conditions(
-        conditions=conditions, start=start, end=end, allow_partial=allow_partial
+        conditions=conditions, start=start, stop=stop, allow_partial=allow_partial
     )
     sql = f"DELETE FROM {table_name}"
     if not where:
@@ -348,7 +403,7 @@ def _select_records_sql(
     conditions: dict,
     columns: typing.Optional[typing.List[str]] = None,
     start: OptionalInt = None,
-    end: OptionalInt = None,
+    stop: OptionalInt = None,
     allow_partial=True,
 ) -> ReturnType:
     """create SQL select statement and values
@@ -361,12 +416,12 @@ def _select_records_sql(
         values to select
     conditions : dict
         the WHERE conditions
-    start, end : OptionalInt
-        select records whose (start, end) values lie between start and end,
+    start, stop : OptionalInt
+        select records whose (start, stop) values lie between start and stop,
         or overlap them if (allow_partial is True)
     allow_partial : bool, optional
-        if False, only records within start, end are included. If True,
-        all records that overlap the segment defined by start, end are included.
+        if False, only records within start, stop are included. If True,
+        all records that overlap the segment defined by start, stop are included.
 
     Returns
     -------
@@ -391,7 +446,7 @@ def _count_records_sql(
     conditions: dict,
     columns: typing.Optional[typing.List[str]] = None,
     start: OptionalInt = None,
-    end: OptionalInt = None,
+    stop: OptionalInt = None,
     allow_partial=True,
 ) -> ReturnType:
     """create SQL count statement and values
@@ -404,12 +459,12 @@ def _count_records_sql(
         values to select
     conditions : dict
         the WHERE conditions
-    start, end : OptionalInt
-        select records whose (start, end) values lie between start and end,
+    start, stop : OptionalInt
+        select records whose (start, stop) values lie between start and stop,
         or overlap them if (allow_partial is True)
     allow_partial : bool, optional
-        if False, only records within start, end are included. If True,
-        all records that overlap the segment defined by start, end are included.
+        if False, only records within start, stop are included. If True,
+        all records that overlap the segment defined by start, stop are included.
 
     Returns
     -------
@@ -448,6 +503,38 @@ def _compatible_schema(db: sqlite3.Connection, schema: dict[str, set]) -> bool:
     return True
 
 
+def _rename_column_if_exists(
+    db: sqlite3.Connection, table_name: str, old_column: str, new_column: str
+) -> None:
+    """Rename a column in a sqlite3 database only if it exists.
+
+    Parameters
+    ----------
+    db : sqlite3.Connection
+        the sqlite3 connection
+    table_name : str
+        the table name
+    old_column : str
+        the column to rename if it exists
+    new_column : str
+        the new name of old_column
+    """
+    cur = db.cursor()
+    table_info = cur.execute(f"PRAGMA table_info({table_name});").fetchall()
+
+    # Check if the column exists
+    for column_info in table_info:
+        if column_info["name"] == old_column:
+            break
+    else:
+        return  # There is no column to rename
+
+    sql = f'ALTER TABLE {table_name} RENAME COLUMN "{old_column}" TO {new_column}'
+
+    cur.execute(sql)
+    db.commit()
+
+
 class SqliteAnnotationDbMixin:
     # table schema for user provided annotations
     _user_schema = {
@@ -456,7 +543,7 @@ class SqliteAnnotationDbMixin:
         "name": "TEXT",
         "parent_id": "TEXT",
         "start": "INTEGER",
-        "end": "INTEGER",
+        "stop": "INTEGER",
         "strand": "TEXT",
         "spans": "array",
         "attributes": "TEXT",
@@ -578,6 +665,9 @@ class SqliteAnnotationDbMixin:
             )
             self._db.row_factory = sqlite3.Row
 
+            for table_name in self.table_names:
+                _rename_column_if_exists(self._db, table_name, "end", "stop")
+
         return self._db
 
     def _execute_sql(self, cmnd, values=None):
@@ -619,10 +709,10 @@ class SqliteAnnotationDbMixin:
             whether the annotation is an alignment annotation
         """
         spans = numpy.array(sorted(sorted(coords) for coords in spans), dtype=int)
-        # the start, end variables are added into the record via the loop
+        # the start, stop variables are added into the record via the loop
         # on local variables
         start = int(spans.min())
-        end = int(spans.max())
+        stop = int(spans.max())
         # we define record as all defined variables from local name space,
         # excluding "self"
         record = {k: v for k, v in locals().items() if k != "self" and v is not None}
@@ -652,7 +742,7 @@ class SqliteAnnotationDbMixin:
         column: str,
         name: str,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
         biotype: OptionalStr = None,
         allow_partial: bool = False,
     ) -> typing.List[FeatureDataType]:
@@ -662,7 +752,7 @@ class SqliteAnnotationDbMixin:
             conditions={column: name, "biotype": biotype},
             columns=columns,
             start=start,
-            end=end,
+            stop=stop,
             allow_partial=allow_partial,
         )
         for result in self._execute_sql(sql, values=vals):
@@ -675,7 +765,7 @@ class SqliteAnnotationDbMixin:
         self, name: str, biotype: OptionalStr = None, **kwargs
     ) -> typing.Iterator[FeatureDataType]:
         """yields children of name"""
-        # kwargs is used because other classes need start / end
+        # kwargs is used because other classes need start / stop
         # just uses search matching
         for table_name in self.table_names:
             columns = "seqid", "biotype", "spans", "strand", "name", "parent_id"
@@ -727,6 +817,11 @@ class SqliteAnnotationDbMixin:
                         )  # remove invalid field for the FeatureDataType
                         yield parent
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_records_matching(
         self,
         *,
@@ -734,7 +829,7 @@ class SqliteAnnotationDbMixin:
         seqid: str = None,
         name: str = None,
         start: int = None,
-        end: int = None,
+        stop: int = None,
         strand: OptionalStr = None,
         attributes: OptionalStr = None,
         on_alignment: bool = None,
@@ -751,6 +846,11 @@ class SqliteAnnotationDbMixin:
             for result in self._get_records_matching(table_name, **kwargs):
                 yield dict(zip(result.keys(), result))
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_features_matching(
         self,
         *,
@@ -758,7 +858,7 @@ class SqliteAnnotationDbMixin:
         seqid: str = None,
         name: str = None,
         start: int = None,
-        end: int = None,
+        stop: int = None,
         strand: OptionalStr = None,
         attributes: OptionalStr = None,
         on_alignment: bool = None,
@@ -1041,6 +1141,11 @@ class SqliteAnnotationDbMixin:
             self.db.backup(backup)
         backup.close()
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def subset(
         self,
         *,
@@ -1049,7 +1154,7 @@ class SqliteAnnotationDbMixin:
         seqid: str = None,
         name: str = None,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
         strand: OptionalStr = None,
         attributes: OptionalStr = None,
         allow_partial: bool = False,
@@ -1057,7 +1162,7 @@ class SqliteAnnotationDbMixin:
         """returns a new db instance with records matching the provided conditions"""
         # make sure python, not numpy, integers
         start = start if start is None else int(start)
-        end = end if end is None else int(end)
+        stop = stop if stop is None else int(stop)
 
         kwargs = {k: v for k, v in locals().items() if k not in {"self", "source"}}
 
@@ -1143,7 +1248,7 @@ class GffAnnotationDb(SqliteAnnotationDbMixin):
         "source": "TEXT",
         "biotype": "TEXT",  # type in GFF
         "start": "INTEGER",
-        "end": "INTEGER",
+        "stop": "INTEGER",
         "score": "TEXT",  # check defn
         "strand": "TEXT",
         "phase": "TEXT",
@@ -1182,12 +1287,12 @@ class GffAnnotationDb(SqliteAnnotationDbMixin):
         # required
         rows = []
         for record in reduced.values():
-            # our Feature code assumes start always < end,
+            # our Feature code assumes start always < stop,
             # we record direction using Strand
             spans = numpy.array(sorted(record["spans"]), dtype=int)  # sorts the rows
             spans.sort(axis=1)
             record["start"] = int(spans.min())
-            record["end"] = int(spans.max())
+            record["stop"] = int(spans.max())
             record["spans"] = spans
             rows.append(tuple(record.get(c) for c in col_order))
 
@@ -1211,6 +1316,7 @@ class GffAnnotationDb(SqliteAnnotationDbMixin):
         # record.
         for record in records:
             record["biotype"] = record.pop("Type")
+            record["stop"] = record.pop("End")
 
             # we force all keys that map to table column names to be lower case
             for key in tuple(record):
@@ -1234,7 +1340,7 @@ class GffAnnotationDb(SqliteAnnotationDbMixin):
                 reduced[record_id]["spans"] = []
 
             # should this just be an append?
-            reduced[record_id]["spans"].append((record["start"], record["end"]))
+            reduced[record_id]["spans"].append((record["start"], record["stop"]))
 
         del records
 
@@ -1267,7 +1373,7 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
         "source": "TEXT",
         "biotype": "TEXT",  # type in GFF
         "start": "INTEGER",
-        "end": "INTEGER",
+        "stop": "INTEGER",
         "strand": "TEXT",
         "comments": "TEXT",
         "spans": "array",  # aggregation of coords across records
@@ -1318,7 +1424,7 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
         rows = []
         exclude = {"translation", "location", "type"}  # location is grabbed directly
         for record in records:
-            # our Feature code assumes start always < end,
+            # our Feature code assumes start always < stop,
             store = {"seqid": seqid}
             # we create the location data directly
             if location := record.get("location", None):
@@ -1326,7 +1432,7 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
                 if strand := location.strand:
                     store["strand"] = "-" if strand == -1 else "+"
                 store["start"] = int(store["spans"].min())
-                store["end"] = int(store["spans"].max())
+                store["stop"] = int(store["spans"].max())
 
             attrs_keys = record.keys() - exclude
             store["biotype"] = record["type"]
@@ -1367,17 +1473,22 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
         self._num_fakeids += 1
         return name
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_feature_children(
         self,
         name: str,
         biotype: OptionalStr = None,
         exclude_biotype: OptionalStr = None,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
     ) -> typing.Iterator[FeatureDataType]:
         """yields children of name"""
         # children must lie within the parent coordinates
-        if start is None or end is None:
+        if start is None or stop is None:
             name = self.__class__.__name__
             msg = f"coordinates required to query children for {name!r}"
             raise ValueError(msg)
@@ -1391,7 +1502,7 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
                 "name",
                 "parent_id",
                 "start",
-                "end",
+                "stop",
             )
             if table_name == "user":
                 columns += ("on_alignment",)
@@ -1403,27 +1514,32 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
                 name=name,
                 biotype=biotype,
                 start=start,
-                end=end,
+                stop=stop,
                 allow_partial=False,
             ):
                 if feat["biotype"] == exclude_biotype:
                     continue
-                cstart, cend = feat.pop("start"), feat.pop("end")
-                if not (start <= cstart < end and start < cend <= end):
+                cstart, cstop = feat.pop("start"), feat.pop("stop")
+                if not (start <= cstart < stop and start < cstop <= stop):
                     continue
 
                 feat.pop("parent_id")  # remove invalid field for the FeatureDataType
                 yield feat
 
+    @c3warn.deprecated_args(
+        "2024.9",
+        reason="replace usage of SQL keyword as column name",
+        old_new=[("end", "stop")],
+    )
     def get_feature_parent(
         self,
         name: str,
         exclude_biotype: OptionalStr = None,
         start: OptionalInt = None,
-        end: OptionalInt = None,
+        stop: OptionalInt = None,
     ) -> typing.Iterator[FeatureDataType]:
         """yields parents of name"""
-        if start is None or end is None:
+        if start is None or stop is None:
             name = self.__class__.__name__
             msg = f"coordinates required to query parent for {name!r}"
             raise ValueError(msg)
@@ -1438,7 +1554,7 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
                 "name",
                 "parent_id",
                 "start",
-                "end",
+                "stop",
             )
             if table_name == "user":
                 columns += ("on_alignment",)
@@ -1448,15 +1564,15 @@ class GenbankAnnotationDb(SqliteAnnotationDbMixin):
                 column="name",
                 name=name,
                 start=start,
-                end=end,
+                stop=stop,
                 allow_partial=False,
             ):
                 # add support for != operation to SQL where clause generation
                 if feat["biotype"] == exclude_biotype:
                     continue
 
-                cstart, cend = feat.pop("start"), feat.pop("end")
-                if cstart > start or end > cend:
+                cstart, cstop = feat.pop("start"), feat.pop("stop")
+                if cstart > start or stop > cstop:
                     continue
 
                 feat.pop("parent_id")  # remove invalid field for the FeatureDataType
@@ -1609,3 +1725,96 @@ def load_annotations(
             show_progress=show_progress,
         )
     )
+
+
+def _update_array_format(data: bytes) -> bytes:
+    """Convert from the old to the new sqlite numpy array format.
+
+    Converts the previous format saved with numpy.ndarray.tobytes
+    to the .npy format generated with numpy.save.
+
+    Ignores any entries saved in the new format
+
+    Parameters
+    ----------
+    data : bytes
+        The sqlite3 representation of the numpy array.
+
+    Returns
+    -------
+    bytes
+        The new sqlite3 representation of the numpy array.
+    """
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        new_data = array_to_sqlite(sqlite_to_array(data))
+    return new_data
+
+
+def update_file_format(
+    source_path: os.PathLike,
+    db_class: typing.Union[
+        type[BasicAnnotationDb],
+        type[GenbankAnnotationDb],
+        type[GffAnnotationDb],
+    ],
+    backup: bool = True,
+) -> None:
+    """Update the database file to the latest format.
+
+    Fixes an OS dependent issue under the previous representation.
+
+    Ensure update_file_format is run on the same OS used to generate the
+    database file.
+
+    By default performs a backup of the database to another file before
+    updating the given file to the new format. This is in case the function
+    is run on a different OS than the one used to generate the file where
+    the spans column may become corrupted. The backup file is saved as the
+    original file path appended with ".bak".
+
+    See https://github.com/cogent3/cogent3/issues/1776 for more details.
+
+    Parameters
+    ----------
+    source_path : os.PathLike
+        The database file to reformat.
+    db_class : typing.Union[ type[BasicAnnotationDb], type[GenbankAnnotationDb], type[GffAnnotationDb], ]
+        The type of database the file is.
+    backup : bool, optional
+        If True (default), performs a backup of the database before updating.
+        Otherwise does not perform a backup prior to update (not recommended).
+    """
+    source_path = pathlib.Path(source_path).expanduser()
+
+    if not source_path.exists():
+        raise OSError(f"File {source_path} does not exist.")
+
+    anno_db = db_class(source=source_path)
+
+    if backup:
+        backup_path = source_path.parent / f"{source_path.name}.bak"
+        if backup_path.exists():
+            raise FileExistsError(
+                f"Backup file already exists for {source_path}. "
+                f"If there was a problem with the conversion process, "
+                f"update_file_format should be run on the backed up file. "
+                f"Ensure update_file_format is run on the same OS used to generate the file."
+            )
+        anno_db.write(backup_path)
+
+    conn = anno_db.db
+    conn.create_function("update_array_format", 1, _update_array_format)
+    for table_name in anno_db.table_names:
+        cursor = conn.cursor()
+        array_columns = [
+            r["name"]
+            for r in cursor.execute(f"PRAGMA table_info({table_name});").fetchall()
+            if r["type"] == "array"
+        ]
+
+        for column in array_columns:
+            cursor.execute(
+                f"UPDATE {table_name} SET {column}=update_array_format({column});"
+            )
+        conn.commit()
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index b2e782ccb..022fd5a57 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -10,6 +10,7 @@ Sequences are intended to be immutable. This is not enforced by the code for
 performance reasons, but don't alter the MolType or the sequence data after
 creation.
 """
+
 from __future__ import annotations
 
 import contextlib
@@ -920,7 +921,7 @@ class Sequence(SequenceI):
         name
             name of the feature
         start, stop
-            start, end positions to search between, relative to offset
+            start, stop positions to search between, relative to offset
             of this sequence. If not provided, entire span of sequence is used.
 
         Notes
@@ -955,11 +956,11 @@ class Sequence(SequenceI):
         query_start = self._seq.absolute_position(start, include_boundary=False)
         # we set include_boundary=True because stop is exclusive indexing,
         # i,e., the stop can be equal to the length of the view
-        query_end = self._seq.absolute_position(stop, include_boundary=True)
+        query_stop = self._seq.absolute_position(stop, include_boundary=True)
 
         rev_strand = strand == -1
         if rev_strand:
-            query_start, query_end = query_end, query_start
+            query_start, query_stop = query_stop, query_start
 
         query_start = max(query_start, 0)
         # in the underlying db, features are always plus strand oriented
@@ -983,7 +984,7 @@ class Sequence(SequenceI):
             name=name,
             biotype=biotype,
             start=query_start,
-            end=query_end,  # todo gah end should be stop
+            stop=query_stop,
             allow_partial=allow_partial,
         ):
             # spans need to be converted from absolute to relative positions
diff --git a/src/cogent3/evolve/likelihood_function.py b/src/cogent3/evolve/likelihood_function.py
index 64aed6862..46225c116 100644
--- a/src/cogent3/evolve/likelihood_function.py
+++ b/src/cogent3/evolve/likelihood_function.py
@@ -1,5 +1,6 @@
 import json
 import random
+import typing
 
 from collections import defaultdict
 from copy import deepcopy
@@ -8,6 +9,7 @@ import numpy
 
 from cogent3._version import __version__
 from cogent3.core.alignment import ArrayAlignment
+from cogent3.core.tree import PhyloNode
 from cogent3.evolve import substitution_model
 from cogent3.evolve.simulate import AlignmentEvolver, random_sequence
 from cogent3.maths.matrix_exponential_integration import expected_number_subs
@@ -17,6 +19,7 @@ from cogent3.maths.measure import (
     paralinear_discrete_time,
 )
 from cogent3.recalculation.definition import ParameterController
+from cogent3.recalculation.scope import InvalidScopeError
 from cogent3.util import table
 from cogent3.util.dict_array import DictArrayTemplate
 from cogent3.util.misc import adjusted_gt_minprob, get_object_provenance
@@ -376,10 +379,8 @@ class LikelihoodFunction(ParameterController):
                 length = self.get_param_value("length", edge=name, **kw)
                 array *= length
         except KeyError as err:
-            if err[0] == "Q" and name != "Q":
-                raise RuntimeError("rate matrix not known by this model")
-            else:
-                raise
+            raise InvalidScopeError from err
+
         return DictArrayTemplate(self._motifs, self._motifs).wrap(array)
 
     def _getLikelihoodValuesSummedAcrossAnyBins(self, locus=None):
@@ -564,7 +565,7 @@ class LikelihoodFunction(ParameterController):
         results = [title, lnL, nfp] + results if lnL else [title, nfp] + results
         return "\n".join(map(str, results))
 
-    def get_annotated_tree(self, length_as=None):
+    def get_annotated_tree(self, length_as: typing.Optional[str] = None) -> PhyloNode:
         """returns tree with model attributes on node.params
 
         length_as : str or None
@@ -582,15 +583,15 @@ class LikelihoodFunction(ParameterController):
 
         is_discrete = isinstance(self.model, DiscreteSubstitutionModel)
 
-        if is_discrete and not length_as == "paralinear":
+        if is_discrete and length_as != "paralinear":
             raise ValueError(f"{length_as} invalid for discrete time process")
 
         assert length_as in ("ENS", "paralinear", None)
         d = self.get_param_value_dict(["edge"])
         lengths = d.pop("length", None)
         mprobs = self.get_motif_probs_by_node()
-        if not is_discrete:
-            ens = self.get_lengths_as_ens(motif_probs=mprobs)
+
+        ens = {} if is_discrete else self.get_lengths_as_ens(motif_probs=mprobs)
 
         plin = self.get_paralinear_metric(motif_probs=mprobs)
         if length_as == "ENS":
@@ -604,20 +605,46 @@ class LikelihoodFunction(ParameterController):
                 edge.params["mprobs"] = mprobs[edge.name].to_dict()
                 continue
 
-            if not is_discrete:
-                edge.params["ENS"] = ens[edge.name]
-
+            edge.params["ENS"] = ens.get(edge.name)
             edge.params["length"] = lengths[edge.name]
             edge.params["paralinear"] = plin[edge.name]
             edge.params["mprobs"] = mprobs[edge.name].to_dict()
             for par in d:
-                val = d[par][edge.name]
+                val = d[par].get(edge.name)
                 if par == length_as:
                     val = ens[edge.name]
                 edge.params[par] = val
 
         return tree
 
+    def get_ens_tree(self) -> PhyloNode:
+        """returns tree with length as ENS
+
+        Notes
+        -----
+        The paralinear distance is added to node.params["paralinear"].
+
+        If it's a discrete-time model, branch lengths are set to None.
+
+        For a stationary model, branch lengths will be unchanged from
+        those values displayed in the statistics tables.
+        """
+        from cogent3.evolve.ns_substitution_model import (
+            DiscreteSubstitutionModel,
+        )
+
+        mprobs = self.get_motif_probs_by_node()
+        if isinstance(self.model, DiscreteSubstitutionModel):
+            raise TypeError("cannot get ENS for discrete-time models")
+
+        ens = self.get_lengths_as_ens(motif_probs=mprobs)
+
+        tree = self._tree.deepcopy()
+        for edge in tree.get_edge_vector(include_root=False):
+            edge.params["length"] = ens[edge.name]
+
+        return tree
+
     def get_motif_probs(self, edge=None, bin=None, locus=None, position=None):
         """
         Parameters
@@ -702,23 +729,22 @@ class LikelihoodFunction(ParameterController):
         motif_probs : dict or DictArray
             an item for each edge of the tree. Computed if not provided.
         """
-        from cogent3.evolve.ns_substitution_model import (
-            DiscreteSubstitutionModel,
-        )
-
-        is_discrete = isinstance(self.model, DiscreteSubstitutionModel)
-
         if motif_probs is None:
             motif_probs = self.get_motif_probs_by_node()
+
         plin = {}
         for edge in self.tree.get_edge_vector(include_root=False):
             parent_name = edge.parent.name
             pi = motif_probs[parent_name]
             P = self.get_psub_for_edge(edge.name)
-            if is_discrete:
+            try:
+                Q = self.get_rate_matrix_for_edge(edge.name, calibrated=False)
+            except InvalidScopeError:
+                Q = None
+
+            if Q is None:
                 para = paralinear_discrete_time(P.array, pi.array)
             else:
-                Q = self.get_rate_matrix_for_edge(edge.name, calibrated=False)
                 para = paralinear_continuous_time(P.array, pi.array, Q.array)
 
             plin[edge.name] = para
@@ -729,16 +755,31 @@ class LikelihoodFunction(ParameterController):
         """returns {edge.name: ens, ...} where ens is the expected number of substitutions
 
         for a stationary Markov process, this is just branch length
+
         Parameters
         ----------
         motif_probs : dict or DictArray
             an item for each edge of the tree. Computed if not provided.
         """
+        from cogent3.evolve.ns_substitution_model import (
+            DiscreteSubstitutionModel,
+        )
+
         if motif_probs is None:
             motif_probs = self.get_motif_probs_by_node()
 
         edge_parent = self.tree.child_parent_map()
-        lengths = {e: self.get_param_value("length", edge=e) for e in edge_parent}
+        lengths = {}
+        for e in edge_parent:
+            try:
+                length = self.get_param_value("length", edge=e)
+            except (InvalidScopeError, KeyError):
+                length = None
+            lengths[e] = length
+
+        if isinstance(self.model, DiscreteSubstitutionModel):
+            return lengths
+
         if not isinstance(self.model, substitution_model.Stationary):
             ens = {}
             for e in edge_parent:
diff --git a/src/cogent3/maths/matrix_exponential_integration.py b/src/cogent3/maths/matrix_exponential_integration.py
index 905adbe10..c87246388 100644
--- a/src/cogent3/maths/matrix_exponential_integration.py
+++ b/src/cogent3/maths/matrix_exponential_integration.py
@@ -58,15 +58,17 @@ class VanLoanIntegratingExponentiator(_Exponentiator):
 
 
 class VonBingIntegratingExponentiator(_Exponentiator):
-
     """An exponentiator that evaluates int_0^t exp(Q*s)ds
-    using the method of Von Bing."""
+    using the method of Von Bing Yap (Personal Communication)."""
 
     def __init__(self, Q):
         """
-        Q -- a diagonisable matrix.
+        Parameters
+        ----------
+        Q
+            a diagonisable matrix.
         """
-        super(VonBingIntegratingExponentiator, self).__init__(Q)
+        super().__init__(Q)
         self.roots, self.evT = eig(Q)
         self.evI = inv(self.evT.T)
         # Remove following check if performance is a concern
diff --git a/src/cogent3/maths/stats/distribution.py b/src/cogent3/maths/stats/distribution.py
index aad172bd9..7a891b69e 100644
--- a/src/cogent3/maths/stats/distribution.py
+++ b/src/cogent3/maths/stats/distribution.py
@@ -23,6 +23,7 @@ from cogent3.maths.stats.special import (
     log1p,
     ndtri,
 )
+from cogent3.util import warning as c3warn
 
 
 # ndtri import b/c it should be available via this module
@@ -43,27 +44,33 @@ def tprob(x, df):
     return 2 * t.sf(abs(x), df)
 
 
-def poisson_high(successes, mean):
-    """Returns right tail of Poission distribution, Pr(X > x).
-
-    successes ranges from 0 to infinity. mean must be positive.
-    """
+@c3warn.deprecated_callable(
+    version="2024.9",
+    reason="use scipy.stats.poisson.sf() instead",
+    is_discontinued=True,
+)
+def poisson_high(successes, mean):  # pragma: no cover
+    """being removed"""
     return pdtrc(successes, mean)
 
 
-def poisson_low(successes, mean):
-    """Returns left tail of Poisson distribution, Pr(X <= x).
-
-    successes ranges from 0 to infinity. mean must be positive.
-    """
+@c3warn.deprecated_callable(
+    version="2024.9",
+    reason="use scipy.stats.poisson.cdf() instead",
+    is_discontinued=True,
+)
+def poisson_low(successes, mean):  # pragma: no cover
+    """being removed"""
     return pdtr(successes, mean)
 
 
-def poisson_exact(successes, mean):
-    """Returns Poisson probablity for exactly Pr(X=successes).
-
-    Formula is e^-(mean) * mean^(successes) / (successes)!
-    """
+@c3warn.deprecated_callable(
+    version="2024.9",
+    reason="use scipy.stats.poisson.pmf() instead",
+    is_discontinued=True,
+)
+def poisson_exact(successes, mean):  # pragma: no cover
+    """being removed"""
     if successes == 0:
         return pdtr(0, mean)
     elif successes < mean:  # use left tail
diff --git a/src/cogent3/maths/stats/special.py b/src/cogent3/maths/stats/special.py
index 18a8b2a33..09a4cf37c 100644
--- a/src/cogent3/maths/stats/special.py
+++ b/src/cogent3/maths/stats/special.py
@@ -1,9 +1,10 @@
 """Translations of functions from Release 2.3 of the Cephes Math Library,
 (c) Stephen L. Moshier 1984, 1995.
 """
-
 from numpy import exp, floor, log, sin, sqrt
 
+from cogent3.util import warning as c3warn
+
 
 log_epsilon = 1e-6  # for threshold in log/exp close to 1
 # For IEEE arithmetic (IBMPC):
@@ -58,7 +59,10 @@ def one_minus_exp(x):
         return 1 - exp(x)
 
 
-def permutations(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="use scipy.special.perm() instead", is_discontinued=True
+)
+def permutations(n, k):  # pragma: no cover
     """Returns the number of ways of choosing k items from n, in order.
 
     Defined as n!/(n-k)!.
@@ -76,7 +80,10 @@ def permutations(n, k):
         return exp(ln_permutations(n, k))
 
 
-def permutations_exact(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="use scipy.special.perm() instead", is_discontinued=True
+)
+def permutations_exact(n, k):  # pragma: no cover
     """Calculates permutations by integer division.
 
     Preferred method for small permutations, but slow on larger ones.
@@ -89,7 +96,10 @@ def permutations_exact(n, k):
     return product
 
 
-def ln_permutations(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="use scipy.special.perm() instead", is_discontinued=True
+)
+def ln_permutations(n, k):  # pragma: no cover
     """Calculates permutations by difference in log of gamma function.
 
     Preferred method for large permutations, but slow on smaller ones.
@@ -99,7 +109,10 @@ def ln_permutations(n, k):
     return lgam(n + 1) - lgam(n - k + 1)
 
 
-def combinations(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="Use scipy.special.comb() instead", is_discontinued=True
+)
+def combinations(n, k):  # pragma: no cover
     """Returns the number of ways of choosing k items from n, in order.
 
     Defined as n!/(k!(n-k)!).
@@ -118,7 +131,10 @@ def combinations(n, k):
         return exp(ln_combinations(n, k))
 
 
-def combinations_exact(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="Use scipy.special.comb() instead", is_discontinued=True
+)
+def combinations_exact(n, k):  # pragma: no cover
     """Calculates combinations by integer division.
 
     Preferred method for small combinations, but slow on larger ones.
@@ -146,7 +162,10 @@ def combinations_exact(n, k):
     return product
 
 
-def ln_combinations(n, k):
+@c3warn.deprecated_callable(
+    version="2024.9", reason="Use scipy.special.comb() instead", is_discontinued=True
+)
+def ln_combinations(n, k):  # pragma: no cover
     """Calculates combinations by difference in log of gamma function.
 
     Preferred method for large combinations, but slow on smaller ones.
@@ -167,12 +186,15 @@ def ln_binomial(successes, trials, prob):
 
     Note: no error checking (expects to be called through binomial_exact())
     """
+    from scipy.special import loggamma
+
     prob = fix_rounding_error(prob)
-    return (
-        ln_combinations(trials, successes)
-        + successes * log(prob)
-        + (trials - successes) * log(1.0 - prob)
+    ln_comb = (
+        loggamma(trials + 1)
+        - loggamma(successes + 1)
+        - loggamma(trials - successes + 1)
     )
+    return ln_comb + successes * log(prob) + (trials - successes) * log(1.0 - prob)
 
 
 # Translations of functions from Cephes Math Library, by Stephen L. Moshier
@@ -248,16 +270,26 @@ ZU = [
 ]
 
 
-def erf(a):
-    """Returns the error function of a: see Cephes docs."""
+@c3warn.deprecated_callable(
+    version="2024.9",
+    reason="use scipy.special.erf(x, out=None) instead",
+    is_discontinued=True,
+)
+def erf(a):  # pragma: no cover
+    """being removed"""
     if abs(a) > 1:
         return 1 - erfc(a)
     z = a * a
     return a * polevl(z, ZT) / polevl(z, ZU)
 
 
-def erfc(a):
-    """Returns the complement of the error function of a: see Cephes docs."""
+@c3warn.deprecated_callable(
+    version="2024.9",
+    reason="use scipy.special.erfc(z, out=None) instead",
+    is_discontinued=True,
+)
+def erfc(a):  # pragma: no cover
+    """being removed"""
     if a < 0:
         x = -a
     else:
diff --git a/src/cogent3/util/warning.py b/src/cogent3/util/warning.py
index 5994a8205..f4859b831 100644
--- a/src/cogent3/util/warning.py
+++ b/src/cogent3/util/warning.py
@@ -203,6 +203,11 @@ def deprecated_callable(
         sig = set(inspect.signature(func).parameters)
         _type = "method" if sig & {"self", "cls", "klass"} else "function"
         old = func.__name__
+        if is_discontinued and old == "__init__":
+            # we're really deprecating a class, so get that name
+            old = func.__qualname__.split(".")[-2]
+            _type = "class"
+
         params = dict(
             _type=_type,
             old=old,


Test patch:
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index 838daaed4..84c79fa00 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -1019,15 +1019,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
         got = aln.filtered(func, motif_length=2, drop_remainder=True, warn=False)
         self.assertEqual(len(got), 4)
 
-    def test_positions(self):
-        """SequenceCollection positions property should iterate over positions, using self.names"""
-        r = self.Class({"a": "AAAAAA", "b": "AAA---", "c": "AAAA--"})
-        r.names = ["a", "b", "c"]
-        self.assertEqual(
-            list(r.positions),
-            list(map(list, ["AAA", "AAA", "AAA", "A-A", "A--", "A--"])),
-        )
-
     def test_take_positions_info(self):
         aln = self.Class(
             {"a": "AAAAAAA", "b": "A--A-AA", "c": "AA-----"}, info={"key": "value"}
@@ -1035,36 +1026,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
         tps = aln.take_positions([5, 4, 0])
         self.assertEqual(tps.info["key"], "value")
 
-    def test_get_position_indices(self):
-        """SequenceCollection get_position_indices should return names of cols where f(col)"""
-
-        def gap_1st(x):
-            return x[0] == "-"
-
-        def gap_2nd(x):
-            return x[1] == "-"
-
-        def gap_3rd(x):
-            return x[2] == "-"
-
-        def is_list(x):
-            return isinstance(x, list)
-
-        self.gaps = self.Class(self.gaps.named_seqs, names=["a", "b", "c"])
-
-        self.assertEqual(self.gaps.get_position_indices(gap_1st), [])
-        self.assertEqual(self.gaps.get_position_indices(gap_2nd), [1, 2, 4])
-        self.assertEqual(self.gaps.get_position_indices(gap_3rd), [2, 3, 4, 5, 6])
-        self.assertEqual(self.gaps.get_position_indices(is_list), [0, 1, 2, 3, 4, 5, 6])
-        # should be able to negate
-        self.assertEqual(
-            self.gaps.get_position_indices(gap_2nd, negate=True), [0, 3, 5, 6]
-        )
-        self.assertEqual(
-            self.gaps.get_position_indices(gap_1st, negate=True), [0, 1, 2, 3, 4, 5, 6]
-        )
-        self.assertEqual(self.gaps.get_position_indices(is_list, negate=True), [])
-
     def test_take_positions_if(self):
         """SequenceCollection take_positions_if should return cols where f(col) is True"""
 
@@ -1093,6 +1054,7 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
             {"a": "AAAAA", "b": "-A-AA", "c": "-----"},
         )
         self.assertEqual(self.gaps.take_positions_if(is_list), self.gaps)
+        self.assertEqual(self.gaps.take_positions_if(is_list), self.gaps)
 
         self.assertTrue(
             isinstance(self.gaps.take_positions_if(gap_1st), _SequenceCollectionBase)
@@ -2025,8 +1987,8 @@ class ArrayAlignmentSpecificTests(TestCase):
     def test_iter(self):
         """ArrayAlignment iter should iterate over positions"""
         result = list(iter(self.a2))
-        for i, j in zip(result, [list(i) for i in ["AD", "BE", "CF"]]):
-            self.assertEqual(i, j)
+        for i, j in zip(result, ["AD", "BE", "CF"]):
+            self.assertEqual(i, list(j))
 
     def test_getitem(self):
         """ArrayAlignment getitem act like standard alignment slice"""
@@ -3612,3 +3574,46 @@ def test_to_json(cls):
     got = json.loads(txt)
     expect = aln.to_rich_dict()
     assert got == expect
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_get_position_indices(cls):
+    """get_position_indices should return names of cols where f(col)"""
+
+    def gap_1st(x):
+        return x[0] == "-"
+
+    def gap_2nd(x):
+        return x[1] == "-"
+
+    def gap_3rd(x):
+        return x[2] == "-"
+
+    def is_list(x):
+        return isinstance(x, list)
+
+    gaps = cls({"a": "AAAAAAA", "b": "A--A-AA", "c": "AA-----"}, names=["a", "b", "c"])
+
+    assert gaps.get_position_indices(gap_1st) == []
+    assert gaps.get_position_indices(gap_2nd) == [1, 2, 4]
+    assert gaps.get_position_indices(gap_3rd) == [2, 3, 4, 5, 6]
+    assert gaps.get_position_indices(is_list) == [0, 1, 2, 3, 4, 5, 6]
+    # should be able to negate
+    assert gaps.get_position_indices(gap_2nd, negate=True) == [0, 3, 5, 6]
+    assert gaps.get_position_indices(gap_1st, negate=True) == [0, 1, 2, 3, 4, 5, 6]
+    assert gaps.get_position_indices(is_list, negate=True) == []
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_positions(cls):
+    """positions property should iterate over positions, using self.names"""
+    r = cls({"a": "AAAAAA", "b": "AAA---", "c": "AAAA--"})
+    r.names = ["a", "b", "c"]
+    expect = [list(v) for v in ("AAA", "AAA", "AAA", "A-A", "A--", "A--")]
+    assert list(r.positions) == expect
+
+
+def test_array_align_error_with_mixed_length():
+    data = dict(s1="ACGG", s2="A-G")
+    with pytest.raises(ValueError, match=".* not all the same length.*"):
+        make_aligned_seqs(data=data)
diff --git a/tests/test_core/test_alphabet.py b/tests/test_core/test_alphabet.py
index a59151dd9..d7168a0ae 100644
--- a/tests/test_core/test_alphabet.py
+++ b/tests/test_core/test_alphabet.py
@@ -6,6 +6,9 @@ import pickle
 
 from unittest import TestCase
 
+import numpy
+import pytest
+
 from numpy import unravel_index
 from numpy.testing import assert_equal
 
@@ -304,3 +307,18 @@ class JointEnumerationTests(TestCase):
         v = [7, 15, 18, 0]
         result = a.unpack_arrays(v)
         assert_equal(result, array([[0, 1, 2, 0], [3, 3, 1, 0], [1, 1, 0, 0]]))
+
+
+@pytest.mark.parametrize("data", ("AUGGA", b"AUGGA"))
+def test_to_indices(data):
+    alphabet = RNA.alphabets.degen_gapped
+    got = alphabet.to_indices(data)
+    assert got.tolist() == [2, 0, 3, 3, 2]
+
+
+@pytest.mark.parametrize("cls", (list, numpy.array, tuple))
+def test_from_indices(cls):
+    data = cls([2, 0, 3, 3, 2])
+    alphabet = RNA.alphabets.degen_gapped
+    got = alphabet.from_indices(data)
+    assert got == "AUGGA"
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index 8a4e815ba..05eb53188 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -1,3 +1,6 @@
+import io
+import warnings
+
 import numpy
 import pytest
 
@@ -8,7 +11,9 @@ from cogent3.core.annotation_db import (
     GffAnnotationDb,
     SupportsFeatures,
     _matching_conditions,
+    _rename_column_if_exists,
     load_annotations,
+    update_file_format,
 )
 from cogent3.core.sequence import Sequence
 from cogent3.parse.genbank import MinimalGenbankParser
@@ -270,7 +275,7 @@ def test_empty_data():
 
 
 # testing GenBank files
-@pytest.fixture(scope="session")
+@pytest.fixture(scope="function")
 def gb_db(DATA_DIR):
     return load_annotations(path=DATA_DIR / "annotated_seq.gb")
 
@@ -292,7 +297,7 @@ def test_gb_get_children(gb_db, parent_biotype, name):
             name=name,
             exclude_biotype=parent_biotype,
             start=coords.min(),
-            end=coords.max(),
+            stop=coords.max(),
         )
     )[0]
     assert child["biotype"] != parent["biotype"]
@@ -308,7 +313,7 @@ def test_gb_get_parent(gb_db):
             name=cds_id,
             exclude_biotype="CDS",
             start=coords.min(),
-            end=coords.max(),
+            stop=coords.max(),
         )
     )[0]
     assert parent["biotype"] != cds["biotype"]
@@ -476,8 +481,8 @@ def test_get_features_matching_start_stop(DATA_DIR, seq):
 
 
 def test_matching_conditions():
-    got, _ = _matching_conditions({"start": 1, "end": 5}, allow_partial=True)
-    expect = "((start >= 1 AND end <= 5) OR (start <= 1 AND end > 1) OR (start < 5 AND end >= 5) OR (start <= 1 AND end >= 5))"
+    got, _ = _matching_conditions({"start": 1, "stop": 5}, allow_partial=True)
+    expect = "((start >= 1 AND stop <= 5) OR (start <= 1 AND stop > 1) OR (start < 5 AND stop >= 5) OR (start <= 1 AND stop >= 5))"
     assert got == expect
 
 
@@ -921,6 +926,165 @@ def test_write(gb_db, tmp_path):
     assert isinstance(got, GenbankAnnotationDb)
 
 
+def convert_to_old_np_format(data: bytes) -> bytes:
+    with io.BytesIO(data) as stream:
+        output = numpy.load(stream).tobytes()
+    return output
+
+
+def convert_spans_column(db, table_name):
+    # Convert the database to the old numpy format.
+    conn = db.db
+    conn.create_function("convert_to_old_np_format", 1, convert_to_old_np_format)
+    cursor = conn.cursor()
+    cursor.execute(f"UPDATE {table_name} SET spans = convert_to_old_np_format(spans);")
+    conn.commit()
+
+
+@pytest.mark.parametrize(
+    "db_name,cls", (("gb_db", GenbankAnnotationDb), ("gff_db", GffAnnotationDb))
+)
+def test_read_old_format(db_name, cls, tmp_path, request):
+    db = request.getfixturevalue(db_name)
+
+    path = tmp_path / f"old_format_{db_name}.db"
+
+    old_tables = db.to_rich_dict()["tables"]
+
+    table_name = db_name.split("_")[0]
+    convert_spans_column(db, table_name)
+
+    db.write(path)
+
+    new_db = cls(source=path)
+    with pytest.warns(UserWarning):
+        new_tables = new_db.to_rich_dict()["tables"]
+    assert new_tables == old_tables
+
+
+@pytest.mark.parametrize(
+    "db_name,cls,backup",
+    (
+        ("gb_db", GenbankAnnotationDb, True),
+        ("gff_db", GffAnnotationDb, True),
+        ("gb_db", GenbankAnnotationDb, False),
+        ("gff_db", GffAnnotationDb, False),
+    ),
+)
+def test_update_file_format(db_name, cls, backup, tmp_path, request):
+    db = request.getfixturevalue(db_name)
+    old_tables = db.to_rich_dict()["tables"]
+
+    path = tmp_path / f"old_format_{db_name}.db"
+
+    # Convert to old numpy format
+    table_name = db_name.split("_")[0]
+    convert_spans_column(db, table_name)
+    db.write(path)
+
+    with open(path, "rb") as f:
+        old_format_bytes = f.read()
+
+    # The file should be updated without warning
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        update_file_format(path, cls, backup=backup)
+
+    backup_path = tmp_path / f"old_format_{db_name}.db.bak"
+    if backup:
+        # Check the backup hasn't been modified
+        with open(backup_path, "rb") as f:
+            assert old_format_bytes == f.read()
+    else:
+        # A backup file should not exist
+        assert not backup_path.exists()
+
+    # The file should have been modified
+    with open(path, "rb") as f:
+        assert old_format_bytes != f.read()
+
+    # Reading the new file (when spans is loaded) should not produce warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+
+        new_db = cls(source=path)
+        new_tables = new_db.to_rich_dict()["tables"]
+
+    # The tables should remain the same
+    assert old_tables == new_tables
+
+
+@pytest.mark.parametrize(
+    "db_name,cls",
+    (("gb_db", GenbankAnnotationDb), ("gff_db", GffAnnotationDb)),
+)
+def test_update_file_format_twice_fails(db_name, cls, tmp_path, request):
+    db = request.getfixturevalue(db_name)
+
+    path = tmp_path / f"old_format_{db_name}.db"
+
+    # Convert to old numpy format
+    table_name = db_name.split("_")[0]
+    convert_spans_column(db, table_name)
+    db.write(path)
+
+    update_file_format(path, cls, backup=True)
+    with pytest.raises(FileExistsError):
+        update_file_format(path, cls, backup=True)
+
+
+@pytest.mark.parametrize(
+    "db_name,cls,backup",
+    (
+        ("gb_db", GenbankAnnotationDb, True),
+        ("gff_db", GffAnnotationDb, True),
+        ("gb_db", GenbankAnnotationDb, False),
+        ("gff_db", GffAnnotationDb, False),
+    ),
+)
+def test_update_file_format_does_not_modify_correct_file(
+    db_name, cls, backup, tmp_path, request
+):
+    db = request.getfixturevalue(db_name)
+
+    path = tmp_path / f"correct_format_{db_name}.db"
+    db.write(path)
+
+    with open(path, "rb") as f:
+        old_bytes = f.read()
+
+    # Updating the file format should not produce warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        update_file_format(path, cls, backup=backup)
+
+    # The file should have remained the same
+    with open(path, "rb") as f:
+        assert old_bytes == f.read()
+
+
+@pytest.mark.parametrize(
+    "db_name,cls,backup",
+    (
+        ("gb_db", GenbankAnnotationDb, True),
+        ("gff_db", GffAnnotationDb, True),
+        ("gb_db", GenbankAnnotationDb, False),
+        ("gff_db", GffAnnotationDb, False),
+    ),
+)
+def test_update_file_format_fake_path(db_name, cls, backup, tmp_path):
+    path = tmp_path / f"non_existent_{db_name}.db"
+    with pytest.raises(OSError):
+        update_file_format(path, cls, backup=backup)
+
+    # The path should not have been created in the process
+    assert not path.exists()
+    if backup:
+        # The backup file path should also not exist
+        backup_path = tmp_path / f"non_existent_{db_name}.db.bak"
+        assert not backup_path.exists()
+
+
 @pytest.fixture(scope="function")
 def tmp_dir(tmp_path_factory):
     return tmp_path_factory.mktemp("annotations")
@@ -939,6 +1103,48 @@ def test_load_anns_with_write(DATA_DIR, tmp_dir):
     assert got_data["tables"] == expect_data["tables"]
 
 
+def test_gff_end_renamed_to_stop(gff_db, tmp_path):
+    bad_col_path = tmp_path / "bad_col.gffdb"
+
+    correct_rich_dict = gff_db.to_rich_dict()
+    del correct_rich_dict["init_args"]
+
+    for table_name in gff_db.table_names:
+        # Convert in memory db stop column to the old end name
+        _rename_column_if_exists(gff_db.db, table_name, "stop", "end")
+
+    gff_db.write(tmp_path / bad_col_path)
+
+    # Test that end column is renamed to stop on construction
+    loaded_gff_db = GffAnnotationDb(source=bad_col_path)
+
+    new_rich_dict = loaded_gff_db.to_rich_dict()
+    del new_rich_dict["init_args"]
+
+    assert new_rich_dict == correct_rich_dict
+
+
+def test_gb_end_renamed_to_stop(gb_db, tmp_path):
+    bad_col_path = tmp_path / "bad_col.gbdb"
+
+    correct_rich_dict = gb_db.to_rich_dict()
+    del correct_rich_dict["init_args"]
+
+    for table_name in gb_db.table_names:
+        # Convert in memory db stop column to the old end name
+        _rename_column_if_exists(gb_db.db, table_name, "stop", "end")
+
+    gb_db.write(tmp_path / bad_col_path)
+
+    # Test that end column is renamed to stop on construction
+    loaded_gb_db = GenbankAnnotationDb(source=bad_col_path)
+
+    new_rich_dict = loaded_gb_db.to_rich_dict()
+    del new_rich_dict["init_args"]
+
+    assert new_rich_dict == correct_rich_dict
+
+
 def test_gbdb_get_children_fails_no_coords(gb_db):
     with pytest.raises(ValueError):
         _ = list(gb_db.get_feature_children(name="CNA00110"))
@@ -957,7 +1163,7 @@ def test_load_annotations_invalid_path():
 @pytest.mark.parametrize("integer", (int, numpy.int64))
 def test_subset_gff3_db(gff_db, integer):
     subset = gff_db.subset(
-        seqid="I", start=integer(40), end=integer(70), allow_partial=True
+        seqid="I", start=integer(40), stop=integer(70), allow_partial=True
     )
     # manual inspection of the original GFF3 file indicates 7 records
     # BUT the CDS records get merged into a single row
@@ -965,7 +1171,7 @@ def test_subset_gff3_db(gff_db, integer):
 
 
 def test_subset_empty_db(gff_db):
-    subset = gff_db.subset(seqid="X", start=40, end=70, allow_partial=True)
+    subset = gff_db.subset(seqid="X", start=40, stop=70, allow_partial=True)
     # no records
     assert not len(subset)
 
@@ -975,7 +1181,7 @@ def test_subset_gff3_db_with_user(gff_db):
         seqid="I", name="gene-01", biotype="gene", spans=[(23, 43)], strand="+"
     )
     gff_db.add_feature(**record)
-    subset = gff_db.subset(seqid="I", start=40, end=70, allow_partial=True)
+    subset = gff_db.subset(seqid="I", start=40, stop=70, allow_partial=True)
     # manual inspection of the original GFF3 file indicates 7 records
     # BUT the CDS records get merged into a single row
     assert len(subset) == 7
@@ -990,7 +1196,7 @@ def test_subset_gb_db(gb_db):
 def test_subset_gff3_db_source(gff_db, tmp_dir):
     outpath = tmp_dir / "subset.gff3db"
     subset = gff_db.subset(
-        seqid="I", start=40, end=70, allow_partial=True, source=outpath
+        seqid="I", start=40, stop=70, allow_partial=True, source=outpath
     )
     subset.db.close()
 
diff --git a/tests/test_evolve/test_likelihood_function.py b/tests/test_evolve/test_likelihood_function.py
index 0a3d8acc3..11a3fa8be 100644
--- a/tests/test_evolve/test_likelihood_function.py
+++ b/tests/test_evolve/test_likelihood_function.py
@@ -2261,8 +2261,8 @@ def test_simulate_alignment3():
     assert re.sub("[ATCG]", "x", simulated.to_dict()["a"]) == "x??xxxxxx?"
 
 
-def test_get_lengths_as_ens_matches_manual_calc():
-    """lengths as ENS for a non-reversible model matches a manual calculation"""
+@pytest.fixture(scope="session")
+def gn_mod():
     mprobs = numpy.array([0.1, 0.2, 0.3, 0.4])
     # unequal lengths
     lf = GN().make_likelihood_function(make_tree("(a:0.1,b:0.1,c:0.1)"))
@@ -2272,10 +2272,15 @@ def test_get_lengths_as_ens_matches_manual_calc():
     for par_name in lf.model.get_param_list():
         lf.set_param_rule(par_name, init=init)
         init += 0.1
+    return lf
+
 
-    len_dict = lf.get_lengths_as_ens()
+def test_get_lengths_as_ens_matches_manual_calc(gn_mod):
+    """lengths as ENS for a non-reversible model matches a manual calculation"""
+    mprobs = numpy.array([0.1, 0.2, 0.3, 0.4])
+    len_dict = gn_mod.get_lengths_as_ens()
     expect = expected_number_subs(
-        mprobs, lf.get_rate_matrix_for_edge("a", calibrated=False), 1.0
+        mprobs, gn_mod.get_rate_matrix_for_edge("a", calibrated=False), 1.0
     )
     assert_allclose(len_dict["a"], expect)
 
@@ -2289,4 +2294,125 @@ def test_expm_zero_lengths(expm, DATA_DIR):
     sm = GTR()
     lf = sm.make_likelihood_function(tree=tree, expm=expm)
     lf.set_alignment(aln)
-    lf.optimise(max_evaluations=5, show_progress=False)
+    lf.optimise(max_evaluations=5, show_progress=False, limit_action="ignore")
+
+
+def test_lengths_as_ens_model_mix(DATA_DIR):
+    aln = load_aligned_seqs((DATA_DIR / "brca1_5.paml"), moltype="dna")
+    names = ["Mouse", "Human", "HowlerMon"]
+    aln = aln.take_seqs(names)
+    tree = make_tree((f"({names[0]}:0.0,{names[1]}:0.0,{names[2]}:0.0)"))
+    sm = GTR()
+    lf = sm.make_likelihood_function(tree=tree, discrete_edges=["Mouse"])
+    lf.set_alignment(aln)
+    lengths = lf.get_lengths_as_ens()
+    assert lengths["Mouse"] is None
+    assert lengths["Human"] == lengths["HowlerMon"] == 1
+
+
+@pytest.mark.parametrize("discrete_time", (False, True)[1:])
+def test_get_annotated_tree_model_mix(DATA_DIR, discrete_time):
+    aln = load_aligned_seqs((DATA_DIR / "brca1_5.paml"), moltype="dna")
+    names = ["Mouse", "Human", "HowlerMon"]
+    aln = aln.take_seqs(names)
+    tree = make_tree((f"({names[0]}:0.0,{names[1]}:0.0,{names[2]}:0.0)"))
+    sm = GTR()
+    lf = sm.make_likelihood_function(
+        tree=tree, discrete_edges=["Mouse"] if discrete_time else None
+    )
+    lf.set_alignment(aln)
+    lf.optimise(max_evaluations=5, show_progress=False, limit_action="ignore")
+    anno_tree = lf.get_annotated_tree(length_as="ENS")
+    assert isinstance(anno_tree, type(tree))
+
+
+@pytest.fixture(scope="session")
+def mixed_model(DATA_DIR):
+    aln = load_aligned_seqs((DATA_DIR / "brca1_5.paml"), moltype="dna")
+    names = ["Mouse", "Human", "HowlerMon"]
+    aln = aln.take_seqs(names)
+    tree = make_tree((f"({names[0]}:0.0,{names[1]}:0.0,{names[2]}:0.0)"))
+    sm = GTR()
+    lf = sm.make_likelihood_function(tree=tree, discrete_edges=["Mouse"])
+    lf.set_alignment(aln)
+    return lf
+
+
+@pytest.fixture(scope="session")
+def bh_model(DATA_DIR):
+    aln = load_aligned_seqs((DATA_DIR / "brca1_5.paml"), moltype="dna")
+    names = ["Mouse", "Human", "HowlerMon"]
+    aln = aln.take_seqs(names)
+    tree = make_tree((f"({names[0]}:0.0,{names[1]}:0.0,{names[2]}:0.0)"))
+    sm = get_model("BH")
+    lf = sm.make_likelihood_function(tree=tree)
+    lf.set_alignment(aln)
+    return lf
+
+
+def test_rate_matrix_for_edge_mixed_discrete(mixed_model):
+    from cogent3.recalculation.scope import InvalidScopeError
+    from cogent3.util.dict_array import DictArray
+
+    with pytest.raises(InvalidScopeError):
+        mixed_model.get_rate_matrix_for_edge("Mouse")
+
+    Q = mixed_model.get_rate_matrix_for_edge("Human")
+
+    assert isinstance(Q, DictArray)
+
+
+def test_get_paralinear_metric_mixed(mixed_model):
+    # should not fail with a mixed discrete, continuous-time model
+    dists = mixed_model.get_paralinear_metric()
+    # making sure have no None's
+    assert None not in set(dists.values())
+
+
+def test_get_lengths_as_ens_mixed(mixed_model):
+    # should not fail with a mixed discrete, continuous-time model
+    lengths = mixed_model.get_lengths_as_ens()
+    # and the length for discrete-time edge Mouse should be None
+    assert lengths["Mouse"] is None
+    # comparing against default value
+    assert lengths["Human"] == lengths["HowlerMon"] == 1
+
+
+def test_get_lengths_as_ens_discrete(bh_model):
+    # should not fail with a discrete-time model
+    lengths = bh_model.get_lengths_as_ens()
+    assert lengths["Human"] is lengths["HowlerMon"] is lengths["Mouse"] is None
+
+
+def test_get_annotated_tree_mixed(mixed_model):
+    # should not fail with a mixed discrete, continuous-time model
+    tree = mixed_model.get_annotated_tree(length_as="ENS")
+    # and the length for discrete-time edge Mouse should be None
+    mouse = tree.get_node_matching_name("Mouse")
+    assert mouse.length is None
+    assert mouse.params["ENS"] is None
+    assert mouse.params["paralinear"] is not None
+
+
+def test_get_ens_tree_mixed(mixed_model):
+    # should not fail with a mixed discrete, continuous-time model
+    tree = mixed_model.get_ens_tree()
+    # and the length for discrete-time edge Mouse should be None
+    mouse = tree.get_node_matching_name("Mouse")
+    assert mouse.length is None
+
+
+def test_get_ens_tree_gn(gn_mod):
+    # should not fail with a mixed discrete, continuous-time model
+    tree = gn_mod.get_ens_tree()
+    ens = gn_mod.get_lengths_as_ens()
+    tree_lengths = {
+        edge.name: edge.length for edge in tree.postorder() if edge.name != "root"
+    }
+    assert ens == tree_lengths
+
+
+def test_get_ens_tree_discrete(bh_model):
+    # should fail with a discrete-time model
+    with pytest.raises(TypeError):
+        bh_model.get_ens_tree()
diff --git a/tests/test_maths/test_stats/test_distribution.py b/tests/test_maths/test_stats/test_distribution.py
index f55901bf6..3dd7c6d1a 100644
--- a/tests/test_maths/test_stats/test_distribution.py
+++ b/tests/test_maths/test_stats/test_distribution.py
@@ -20,7 +20,6 @@ from cogent3.maths.stats.distribution import (
     pdtr,
     pdtrc,
     pdtri,
-    poisson_exact,
     poisson_high,
     poisson_low,
     probability_points,
@@ -126,70 +125,6 @@ class DistributionsTests(TestCase):
             for x, p in zip(self.values, probs[df]):
                 assert_almost_equal(tprob(x, df), p, decimal=4)
 
-    def test_poisson_low(self):
-        """Lower tail of poisson should match R for integer successes"""
-        # WARNING: Results only guaranteed for integer successes: floating
-        # point _should_ yield reasonable values, but R rounds to int.
-        expected = {
-            (0, 0): 1,
-            (0, 0.75): 0.4723666,
-            (0, 1): 0.3678794,
-            (0, 5): 0.006737947,
-            (0, 113.7): 4.175586e-50,
-            (2, 0): 1,
-            (2, 3): 0.4231901,
-            (2, 17.8): 3.296636e-06,
-            (17, 29.6): 0.008753318,
-            (180, 0): 1,
-            (180, 137.4): 0.999784,
-            (180, 318): 2.436995e-17,
-            (180, 1024): 8.266457e-233,
-        }
-        for key, value in list(expected.items()):
-            assert_allclose(poisson_low(*key), value, rtol=1e-6)
-
-    def test_poisson_high(self):
-        """Upper tail of poisson should match R for integer successes"""
-        # WARNING: Results only guaranteed for integer successes: floating
-        # point _should_ yield reasonable values, but R rounds to int.
-        expected = {
-            (0, 0): 0,
-            (0, 0.75): 0.5276334,
-            (0, 1): 0.6321206,
-            (0, 5): 0.993262,
-            (0, 113.7): 1,
-            (2, 0): 0,
-            (2, 3): 0.5768099,
-            (2, 17.8): 0.9999967,
-            (17, 29.6): 0.9912467,
-            (180, 0): 0,
-            (180, 137.4): 0.0002159856,
-            (180, 318): 1,
-            (180, 1024): 1,
-        }
-        for key, value in list(expected.items()):
-            assert_allclose(poisson_high(*key), value)
-
-    def test_poisson_exact(self):
-        """Poisson exact should match expected values from R"""
-        expected = {
-            (0, 0): 1,
-            (0, 0.75): 0.4723666,
-            (0, 1): 0.3678794,
-            (0, 5): 0.006737947,
-            (0, 113.7): 4.175586e-50,
-            (2, 0): 0,
-            (2, 3): 0.2240418,
-            (2, 17.8): 2.946919e-06,
-            (17, 29.6): 0.004034353,
-            (180, 0): 0,
-            (180, 137.4): 7.287501e-05,
-            (180, 318): 1.067247e-17,
-            (180, 1024): 6.815085e-233,
-        }
-        for key, value in list(expected.items()):
-            assert_allclose(poisson_exact(*key), value, rtol=1e-6)
-
     def test_binomial_series(self):
         """binomial_exact should match values from R on a whole series"""
         expected = list(
diff --git a/tests/test_maths/test_stats/test_special.py b/tests/test_maths/test_stats/test_special.py
index 5056c37fd..22fe28f40 100644
--- a/tests/test_maths/test_stats/test_special.py
+++ b/tests/test_maths/test_stats/test_special.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Unit tests for special functions used in statistics.
 """
 import math
@@ -8,159 +7,19 @@ from unittest import TestCase
 from numpy.testing import assert_allclose
 
 from cogent3.maths.stats.special import (
-    combinations,
-    combinations_exact,
     igami,
     incbi,
     ln_binomial,
-    ln_combinations,
-    ln_permutations,
     log1p,
     log_one_minus,
     ndtri,
     one_minus_exp,
-    permutations,
-    permutations_exact,
 )
 
 
 class SpecialTests(TestCase):
     """Tests miscellaneous functions."""
 
-    def test_permutations(self):
-        """permutations should return expected results"""
-        self.assertEqual(permutations(1, 1), 1)
-        self.assertEqual(permutations(2, 1), 2)
-        self.assertEqual(permutations(3, 1), 3)
-        self.assertEqual(permutations(4, 1), 4)
-        self.assertEqual(permutations(4, 2), 12)
-        self.assertEqual(permutations(4, 3), 24)
-        self.assertEqual(permutations(4, 4), 24)
-        assert_allclose(permutations(300, 100), 3.8807387193009318e239)
-
-    def test_permutations_errors(self):
-        """permutations should raise errors on invalid input"""
-        self.assertRaises(IndexError, permutations, 10, 50)
-        self.assertRaises(IndexError, permutations, -1, 50)
-        self.assertRaises(IndexError, permutations, 10, -5)
-
-    def test_permutations_float(self):
-        """permutations should use gamma function when floats as input"""
-        assert_allclose(permutations(1.0, 1), 1)
-        assert_allclose(permutations(2, 1.0), 2)
-        assert_allclose(permutations(3.0, 1.0), 3)
-        assert_allclose(permutations(4.0, 1), 4)
-        assert_allclose(permutations(4.0, 2.0), 12)
-        assert_allclose(permutations(4.0, 3.0), 24)
-        assert_allclose(permutations(4, 4.0), 24)
-        assert_allclose(permutations(300, 100), 3.8807387193009318e239)
-
-    def test_permutations_range(self):
-        """permutations should increase gradually with increasing k"""
-        start = 5  # permuations(10,5) = 30240
-        end = 6  # permutations(10,6) = 151200
-        step = 0.1
-        lower_lim = 30240
-        upper_lim = 151200
-        previous_value = 30239.9999
-        while start <= end:
-            obs = permutations(10, start)
-            assert lower_lim <= obs <= upper_lim
-            assert obs > previous_value
-            previous_value = obs
-            start += step
-
-    def test_permutations_exact(self):
-        """permutations_exact should return expected results"""
-        assert_allclose(permutations_exact(1, 1), 1)
-        assert_allclose(permutations_exact(2, 1), 2)
-        assert_allclose(permutations_exact(3, 1), 3)
-        assert_allclose(permutations_exact(4, 1), 4)
-        assert_allclose(permutations_exact(4, 2), 12)
-        assert_allclose(permutations_exact(4, 3), 24)
-        assert_allclose(permutations_exact(4, 4), 24)
-        assert_allclose(permutations_exact(300, 100) / 3.8807387193009318e239, 1.0)
-
-    def test_ln_permutations(self):
-        """ln_permutations should return expected results"""
-        assert_allclose(ln_permutations(1, 1), math.log(1))
-        assert_allclose(ln_permutations(2, 1), math.log(2))
-        assert_allclose(ln_permutations(3, 1.0), math.log(3))
-        assert_allclose(ln_permutations(4, 1), math.log(4))
-        assert_allclose(ln_permutations(4.0, 2), math.log(12))
-        assert_allclose(ln_permutations(4, 3.0), math.log(24))
-        assert_allclose(ln_permutations(4, 4), math.log(24))
-        assert_allclose(ln_permutations(300.0, 100), math.log(3.8807387193009318e239))
-
-    def test_combinations(self):
-        """combinations should return expected results when int as input"""
-        self.assertEqual(combinations(1, 1), 1)
-        self.assertEqual(combinations(2, 1), 2)
-        self.assertEqual(combinations(3, 1), 3)
-        self.assertEqual(combinations(4, 1), 4)
-        self.assertEqual(combinations(4, 2), 6)
-        self.assertEqual(combinations(4, 3), 4)
-        self.assertEqual(combinations(4, 4), 1)
-        self.assertEqual(combinations(20, 4), 19 * 17 * 15)
-        assert_allclose(combinations(300, 100), 4.1582514632578812e81)
-
-    def test_combinations_errors(self):
-        """combinations should raise errors on invalid input"""
-        self.assertRaises(IndexError, combinations, 10, 50)
-        self.assertRaises(IndexError, combinations, -1, 50)
-        self.assertRaises(IndexError, combinations, 10, -5)
-
-    def test_combinations_float(self):
-        """combinations should use gamma function when floats as input"""
-        assert_allclose(combinations(1.0, 1.0), 1)
-        assert_allclose(combinations(2.0, 1.0), 2)
-        assert_allclose(combinations(3.0, 1.0), 3)
-        assert_allclose(combinations(4.0, 1.0), 4)
-        assert_allclose(combinations(4.0, 2), 6)
-        assert_allclose(combinations(4, 3.0), 4)
-        assert_allclose(combinations(4.0, 4.0), 1)
-        assert_allclose(combinations(20.0, 4.0), 19 * 17 * 15)
-        assert_allclose(combinations(300, 100.0), 4.1582514632578812e81)
-
-    def test_combinations_range(self):
-        """combinations should decrease gradually with increasing k"""
-        start = 5  # combinations(10,5) = 252
-        end = 6  # combinations(10,6) = 210
-        step = 0.1
-        lower_lim = 210
-        upper_lim = 252
-        previous_value = 252.00001
-        while start <= end:
-            obs = combinations(10, start)
-            assert lower_lim <= obs <= upper_lim
-            assert obs < previous_value
-            previous_value = obs
-            start += step
-
-    def test_combinations_exact(self):
-        """combinations_exact should return expected results"""
-        self.assertEqual(combinations_exact(1, 1), 1)
-        self.assertEqual(combinations_exact(2, 1), 2)
-        self.assertEqual(combinations_exact(3, 1), 3)
-        self.assertEqual(combinations_exact(4, 1), 4)
-        self.assertEqual(combinations_exact(4, 2), 6)
-        self.assertEqual(combinations_exact(4, 3), 4)
-        self.assertEqual(combinations_exact(4, 4), 1)
-        self.assertEqual(combinations_exact(20, 4), 19 * 17 * 15)
-        assert_allclose(combinations_exact(300, 100), 4.1582514632578812e81)
-
-    def test_ln_combinations(self):
-        """ln_combinations should return expected results"""
-        assert_allclose(ln_combinations(1, 1), math.log(1))
-        assert_allclose(ln_combinations(2, 1), math.log(2))
-        assert_allclose(ln_combinations(3, 1), math.log(3))
-        assert_allclose(ln_combinations(4.0, 1), math.log(4))
-        assert_allclose(ln_combinations(4, 2.0), math.log(6))
-        assert_allclose(ln_combinations(4, 3), math.log(4))
-        assert_allclose(ln_combinations(4, 4.0), math.log(1))
-        assert_allclose(ln_combinations(20, 4), math.log(19 * 17 * 15))
-        assert_allclose(ln_combinations(300, 100), math.log(4.1582514632578812e81))
-
     def test_ln_binomial_integer(self):
         """ln_binomial should match R results for integer values"""
         expected = {
diff --git a/tests/test_util/test_warning.py b/tests/test_util/test_warning.py
index e4dbe916d..9363fca25 100644
--- a/tests/test_util/test_warning.py
+++ b/tests/test_util/test_warning.py
@@ -227,3 +227,14 @@ def test_function_deprecated_args_deprecated_callable_chained_decorators(recwarn
     assert any("argument x which will be removed" in warning for warning in warnings)
     assert any("argument b is discontinued" in warning for warning in warnings)
     assert any("function changed is discontinued" in warning for warning in warnings)
+
+
+def test_class_deprecated(recwarn):
+    class fooclass:
+        @deprecated_callable("2023.6", "Improved change function", is_discontinued=True)
+        def __init__(self):
+            ...
+
+    fooclass()
+    warnings = [warning.message.args[0] for warning in recwarn.list]
+    assert any("fooclass is discontinued" in warning for warning in warnings)