_pyrodigal.pyx

# coding: utf-8
# cython: language_level=3, linetrace=True

"""Bindings to Prodigal, an ORF finder for genomes and metagenomes.

Example:
    Pyrodigal can work on any DNA sequence stored in either a text or a
    byte array. To load a sequence from one of the common sequence formats,
    you can use an external dedicated library such as
    `Biopython <https://github.com/biopython/biopython>`_::

        >>> import gzip
        >>> import Bio.SeqIO
        >>> with gzip.open("pyrodigal/tests/data/KK037166.fna.gz", "rt") as f:
        ...     record = Bio.SeqIO.read(f, "fasta")

    Then use Pyrodigal to find the genes in *metagenomic* mode (without
    training first), and then build a map of codon frequencies for each
    gene::

        >>> from collections import Counter
        >>> import pyrodigal
        >>> p = pyrodigal.OrfFinder(meta=True)
        >>> for gene in p.find_genes(record.seq.encode()):
        ...     gene_seq = gene.sequence()
        ...     codon_counter = Counter()
        ...     for i in range(len(gene_seq), 3):
        ...         codon_counter[gene_seq[i:i+3]] += 1
        ...     codon_frequencies = {
        ...         codon:count/(len(gene_seq)//3)
        ...         for codon, count in codon_counter.items()
        ...     }

Caution:
    In Pyrodigal, sequences are assumed to contain only the usual
    nucleotides (A/T/G/C) as lowercase or uppercase letters; any other
    symbol will be treated as an unknown nucleotide. Be careful to remove
    the gap characters if loading sequences from a multiple alignment file.

References:
    - Hyatt D, Chen GL, Locascio PF, Land ML, Larimer FW, Hauser LJ.
      *Prodigal: prokaryotic gene recognition and translation initiation
      site identification.* BMC Bioinformatics. 2010 Mar 8;11:119.
      doi:10.1186/1471-2105-11-119. PMID:20211023. PMCID:PMC2848648.

"""

# ----------------------------------------------------------------------------

from cpython.buffer cimport PyBUF_READ, PyBUF_WRITE
from cpython.bytes cimport PyBytes_FromStringAndSize, PyBytes_AsString
from cpython.exc cimport PyErr_CheckSignals
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
from cpython.memoryview cimport PyMemoryView_FromMemory
from cpython.ref cimport Py_INCREF
from cpython.tuple cimport PyTuple_New, PyTuple_SET_ITEM
from libc.math cimport sqrt, log, pow, fmax, fmin
from libc.stdint cimport int8_t, uint8_t, uintptr_t
from libc.stdio cimport printf
from libc.stdlib cimport abs, malloc, calloc, free, qsort
from libc.string cimport memcpy, memchr, memset, strstr

from pyrodigal.prodigal cimport bitmap, dprog, gene, node, sequence
from pyrodigal.prodigal.metagenomic cimport NUM_META, _metagenomic_bin, initialize_metagenomic_bins
from pyrodigal.prodigal.node cimport _motif, _node, MIN_EDGE_GENE, MIN_GENE, MAX_SAM_OVLP, cross_mask, compare_nodes, stopcmp_nodes
from pyrodigal.prodigal.sequence cimport _mask, node_type, rcom_seq
from pyrodigal.prodigal.training cimport _training
from pyrodigal._unicode cimport *
from pyrodigal._sequence cimport (
    nucleotide,
    _is_a,
    _is_g,
    _is_gc,
    _is_stop,
    _is_start,
    _is_atg,
    _is_ttg,
    _is_gtg,
    _mer_ndx,
    _letters,
    _complement
)
from pyrodigal.impl.generic cimport skippable_generic

IF TARGET_CPU == "x86":
    from pyrodigal.cpu_features.x86 cimport GetX86Info, X86Info
    IF SSE2_BUILD_SUPPORT:
        from pyrodigal.impl.sse cimport skippable_sse
    IF AVX2_BUILD_SUPPORT:
        from pyrodigal.impl.avx cimport skippable_avx
ELIF TARGET_CPU == "arm" or TARGET_CPU == "aarch64":
    IF TARGET_CPU == "arm":
        from pyrodigal.cpu_features.arm cimport GetArmInfo, ArmInfo
    IF NEON_BUILD_SUPPORT:
        from pyrodigal.impl.neon cimport skippable_neon

IF SYS_IMPLEMENTATION_NAME == "pypy":
    cdef int MVIEW_READ  = PyBUF_READ | PyBUF_WRITE
    cdef int MVIEW_WRITE = PyBUF_READ | PyBUF_WRITE
ELSE:
    cdef int MVIEW_READ  = PyBUF_READ
    cdef int MVIEW_WRITE = PyBUF_WRITE

# ----------------------------------------------------------------------------

import warnings
import textwrap
import threading

include "_version.py"

# --- Module-level constants -------------------------------------------------

cdef int    IDEAL_SINGLE_GENOME = 100000
cdef int    MIN_SINGLE_GENOME   = 20000
cdef int    WINDOW              = 120
cdef size_t MIN_MASKS_ALLOC     = 8
cdef size_t MIN_GENES_ALLOC     = 8
cdef size_t MIN_NODES_ALLOC     = 8 * MIN_GENES_ALLOC
cdef set    TRANSLATION_TABLES  = set(range(1, 7)) | set(range(9, 17)) | set(range(21, 26))

_TRANSLATION_TABLES = TRANSLATION_TABLES


# --- Sequence mask ----------------------------------------------------------

cdef class Mask:
    """The coordinates of a masked region.
    """

    @property
    def begin(self):
        return self.mask.begin

    @property
    def end(self):
        return self.mask.end

cdef class Masks:
    """A list of masked regions within a `~pyrodigal.Sequence`.
    """

    def __cinit__(self):
        self.masks = NULL
        self.capacity = 0
        self.length = 0

    def __init__(self):
        self._clear()

    def __dealloc__(self):
        PyMem_Free(self.masks)

    def __copy__(self):
        return self.copy()

    def __len__(self):
        return self.length

    def __getitem__(self, ssize_t index):
        cdef Mask mask
        if index < 0:
            index += <ssize_t> self.length
        if index >= <ssize_t> self.length or index < 0:
            raise IndexError("list index out of range")
        mask = Mask.__new__(Mask)
        mask.owner = self
        mask.mask = &self.masks[index]
        return mask

    def __sizeof__(self):
        return self.capacity * sizeof(_mask) + sizeof(self)

    cdef inline _mask* _add_mask(
        self,
        const int  begin,
        const int  end,
    ) nogil except NULL:
        """Add a single node to the vector, and return a pointer to that node.
        """

        cdef size_t old_capacity = self.capacity
        cdef _mask* mask

        if self.length >= self.capacity:
            self.capacity = MIN_MASKS_ALLOC if self.capacity == 0 else self.capacity*2
            with gil:
                self.masks = <_mask*> PyMem_Realloc(self.masks, self.capacity * sizeof(_mask))
                if self.masks == NULL:
                    raise MemoryError("Failed to reallocate mask array")
            memset(&self.masks[old_capacity], 0, (self.capacity - old_capacity) * sizeof(_mask))

        self.length += 1
        mask = &self.masks[self.length - 1]
        mask.begin = begin
        mask.end = end
        return mask

    cdef int _clear(self) nogil except 1:
        """Remove all masks from the vector.
        """
        cdef size_t old_length
        old_length, self.length = self.length, 0
        memset(self.masks, 0, old_length * sizeof(_mask))

    def clear(self):
        """Remove all masks from the vector.
        """
        with nogil:
            self._clear()

    cpdef Masks copy(self):
        cdef Masks new = Masks.__new__(Masks)
        new.capacity = self.capacity
        new.length = self.length
        new.masks = <_mask*> PyMem_Malloc(new.capacity * sizeof(_mask))
        if new.masks == NULL:
            raise MemoryError("Failed to allocate masks array")
        memcpy(new.masks, self.masks, new.capacity * sizeof(_mask))
        return new


# --- Input sequence ---------------------------------------------------------

cdef class Sequence:
    """A digitized input sequence.

    Attributes:
        gc (`float`): The GC content of the sequence, as a fraction.
        masks (`~pyrodigal.Masks`): A list of masked regions within the
            sequence.

    """

    # --- Class methods ------------------------------------------------------

    @classmethod
    def from_bytes(cls, const unsigned char[:] sequence, bint mask = False):
        """from_bytes(cls, sequence)\n--

        Create a new `Sequence` object from an ASCII-encoded sequence.

        Arguments:
            sequence (`bytes`): The ASCII-encoded sequence to use. Any
                object implementing the *buffer protocol* is supported.
            mask (`bool`): Enable region-masking for spans of unknown
                characters, preventing genes from being built across them.

        """
        cdef int           i
        cdef int           j
        cdef unsigned char letter
        cdef Sequence      seq
        cdef int           gc_count   = 0
        cdef int           mask_begin = -1

        seq = Sequence.__new__(Sequence)
        seq._allocate(sequence.shape[0])
        seq.masks = Masks.__new__(Masks)

        with nogil:
            for i in range(seq.slen):
                letter = sequence[i]
                if letter == b'A' or letter == b'a':
                    seq.digits[i] = nucleotide.A
                elif letter == b'T' or letter == b't':
                    seq.digits[i] = nucleotide.T
                elif letter == b'G' or letter == b'g':
                    seq.digits[i] = nucleotide.G
                    gc_count += 1
                elif letter == b'C' or letter == b'c':
                    seq.digits[i] = nucleotide.C
                    gc_count += 1
                else:
                    seq.digits[i] = nucleotide.N

            if seq.slen > 0:
                seq.gc = (<double> gc_count) / (<double> seq.slen)

            if mask:
                for i in range(seq.slen):
                    if seq.digits[i] == nucleotide.N:
                        if mask_begin == -1:
                            mask_begin = i
                    else:
                        if mask_begin != -1:
                            seq.masks._add_mask(mask_begin, i-1)
                            mask_begin = -1

        return seq

    @classmethod
    def from_string(cls, str sequence, bint mask = False):
        """from_string(cls, sequence)\n--

        Create a new `Sequence` object from a Unicode sequence.

        Arguments:
            sequence (`str`): The Unicode sequence to use.
            mask (`bool`): Enable region-masking for spans of unknown
                characters, preventing genes from being built across them.

        """
        cdef int      i
        cdef int      j
        cdef Py_UCS4  letter
        cdef Sequence seq
        cdef int      kind
        cdef void*    data
        cdef int      gc_count   = 0
        cdef int      mask_begin = -1

        # make sure the unicode string is in canonical form,
        # --> won't be needed anymore in Python 3.12
        IF SYS_VERSION_INFO_MAJOR <= 3 and SYS_VERSION_INFO_MINOR < 12:
            PyUnicode_READY(sequence)

        seq = Sequence.__new__(Sequence)
        seq._allocate(PyUnicode_GET_LENGTH(sequence))
        seq.masks = Masks.__new__(Masks)

        kind = PyUnicode_KIND(sequence)
        data = PyUnicode_DATA(sequence)

        with nogil:
            for i, j in enumerate(range(0, seq.slen * 2, 2)):
                letter = PyUnicode_READ(kind, data, i)
                if letter == u'A' or letter == u'a':
                    seq.digits[i] = nucleotide.A
                elif letter == u'T' or letter == u't':
                    seq.digits[i] = nucleotide.T
                elif letter == u'G' or letter == u'g':
                    seq.digits[i] = nucleotide.G
                    gc_count += 1
                elif letter == u'C' or letter == u'c':
                    seq.digits[i] = nucleotide.C
                    gc_count += 1
                else:
                    seq.digits[i] = nucleotide.N

            if seq.slen > 0:
                seq.gc = (<double> gc_count) / (<double> seq.slen)

            if mask:
                for i in range(seq.slen):
                    if seq.digits[i] == nucleotide.N:
                        if mask_begin == -1:
                            mask_begin = i
                    else:
                        if mask_begin != -1:
                            seq.masks._add_mask(mask_begin, i-1)
                            mask_begin = -1

        return seq

    # --- Magic methods ------------------------------------------------------

    def __cinit__(self):
        self.slen = 0
        self.gc = 0.0
        self.digits = NULL
        self.masks = None

    def __dealloc__(self):
        PyMem_Free(self.digits)

    def __len__(self):
        return self.slen

    def __sizeof__(self):
        return self.slen * sizeof(uint8_t) + sizeof(self)

    def __str__(self):
        cdef int     i
        cdef Py_UCS4 nuc

        IF SYS_VERSION_INFO_MAJOR <= 3 and SYS_VERSION_INFO_MINOR < 7 and SYS_IMPLEMENTATION_NAME == "pypy":
            cdef bytes dna
            cdef void* data
            # create an empty byte buffer that we can write to
            dna = PyBytes_FromStringAndSize(NULL, self.slen)
            data = <void*> PyBytes_AsString(dna)
        ELSE:
            cdef unicode dna
            cdef int     kind
            cdef void*   data
            # create an empty string that we can write to
            dna  = PyUnicode_New(self.slen, 0x7F)
            kind = PyUnicode_KIND(dna)
            data = PyUnicode_DATA(dna)

        with nogil:
            for i in range(self.slen):
                nuc = _letters[self.digits[i]]
                IF SYS_VERSION_INFO_MAJOR <= 3 and SYS_VERSION_INFO_MINOR < 7 and SYS_IMPLEMENTATION_NAME == "pypy":
                    (<char*> data)[i] = nuc
                ELSE:
                    PyUnicode_WRITE(kind, data, i, nuc)

        IF SYS_VERSION_INFO_MAJOR <= 3 and SYS_VERSION_INFO_MINOR < 7 and SYS_IMPLEMENTATION_NAME == "pypy":
            return dna.decode("ascii")
        ELSE:
            return dna

    # --- C interface -------------------------------------------------------

    cdef int _allocate(self, int slen) except 1:
        self.slen = slen
        self.digits = <uint8_t*> PyMem_Malloc(slen * sizeof(uint8_t))
        if self.digits == NULL:
            raise MemoryError()
        with nogil:
            memset(self.digits, 0, slen * sizeof(uint8_t))
        return 0

    cdef char _amino(
        self,
        int i,
        int tt,
        int strand = 1,
        bint is_init = False,
        char unknown_residue = b"X"
    ) nogil:
        cdef uint8_t x0
        cdef uint8_t x1
        cdef uint8_t x2

        if _is_stop(self.digits, self.slen, i, tt, strand):
            return b"*"
        if _is_start(self.digits, self.slen, i, tt, strand) and is_init:
            return b"M"

        if strand == 1:
            x0 = self.digits[i]
            x1 = self.digits[i+1]
            x2 = self.digits[i+2]
        else:
            x0 = _complement[self.digits[self.slen - 1 - i]]
            x1 = _complement[self.digits[self.slen - 2 - i]]
            x2 = _complement[self.digits[self.slen - 3 - i]]

        if x0 == nucleotide.T and x1 == nucleotide.T and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"F"
        if x0 == nucleotide.T and x1 == nucleotide.T and (x2 == nucleotide.A or x2 == nucleotide.G):
            return b"L"
        if x0 == nucleotide.T and x1 == nucleotide.C and x2 != nucleotide.N:
            return b"S"
        if x0 == nucleotide.T and x1 == nucleotide.A and x2 == nucleotide.T:
            return b"Y"
        if x0 == nucleotide.T and x1 == nucleotide.A and x2 == nucleotide.C:
            return b"Y"
        if x0 == nucleotide.T and x1 == nucleotide.A and x2 == nucleotide.A:
            if tt == 6:
                return b"Q"
            elif tt == 14:
                return b"Y"
        if x0 == nucleotide.T and x1 == nucleotide.A and x2 == nucleotide.G:
            if tt == 6 or tt == 15:
                return b"Q"
            elif tt == 22:
                return b"L"
        if x0 == nucleotide.T and x1 == nucleotide.G and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"C"
        if x0 == nucleotide.T and x1 == nucleotide.G and x2 == nucleotide.A:
            return b"G" if tt == 25 else b"W"
        if x0 == nucleotide.T and x1 == nucleotide.G and x2 == nucleotide.G:
            return b"W"
        if x0 == nucleotide.C and x1 == nucleotide.T and (x2 == nucleotide.T or x2 == nucleotide.C or x2 == nucleotide.A):
            return b"T" if tt == 3 else b"L"
        if x0 == nucleotide.C and x1 == nucleotide.T and x2 == nucleotide.G:
            return b"T" if tt == 3 else b"S" if tt == 12 else b"L"
        if x0 == nucleotide.C and x1 == nucleotide.C and x2 != nucleotide.N:
            return b"P"
        if x0 == nucleotide.C and x1 == nucleotide.A and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"H"
        if x0 == nucleotide.C and x1 == nucleotide.A and (x2 == nucleotide.A or x2 == nucleotide.G):
            return b"Q"
        if x0 == nucleotide.C and x1 == nucleotide.G and x2 != nucleotide.N:
            return b"R"
        if x0 == nucleotide.A and x1 == nucleotide.T and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"I"
        if x0 == nucleotide.A and x1 == nucleotide.T and x2 == nucleotide.A:
            return b"M" if tt == 2 or tt == 3 or tt == 5 or tt == 13 or tt == 22 else b"I"
        if x0 == nucleotide.A and x1 == nucleotide.T and x2 == nucleotide.G:
            return b"M"
        if x0 == nucleotide.A and x1 == nucleotide.C and x2 != nucleotide.N:
            return b"T"
        if x0 == nucleotide.A and x1 == nucleotide.A and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"N"
        if x0 == nucleotide.A and x1 == nucleotide.A and x2 == nucleotide.A:
            return b"N" if tt == 9 or tt == 14 or tt == 21 else b"K"
        if x0 == nucleotide.A and x1 == nucleotide.A and x2 == nucleotide.G:
            return b"K"
        if x0 == nucleotide.A and x1 == nucleotide.G and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"S"
        if x0 == nucleotide.A and x1 == nucleotide.G and (x2 == nucleotide.A or x2 == nucleotide.G):
            return b"G" if tt == 13 else b"S" if tt == 5 or tt == 9 or tt == 14 or tt == 21 else b"R"
        if x0 == nucleotide.G and x1 == nucleotide.T and x2 != nucleotide.N:
            return b"V"
        if x0 == nucleotide.G and x1 == nucleotide.C and x2 != nucleotide.N:
            return b"A"
        if x0 == nucleotide.G and x1 == nucleotide.A and (x2 == nucleotide.T or x2 == nucleotide.C):
            return b"D"
        if x0 == nucleotide.G and x1 == nucleotide.A and (x2 == nucleotide.A or x2 == nucleotide.G):
            return b"E"
        if x0 == nucleotide.G and x1 == nucleotide.G and x2 != nucleotide.N:
            return b"G"

        return unknown_residue

    cdef int _shine_dalgarno_exact(
        self,
        const int pos,
        const int start,
        const _training* tinf,
        const int strand
    ) nogil:
        cdef int i
        cdef int j
        cdef int k
        cdef int mism
        cdef int rdis
        cdef int limit
        cdef int max_val
        cdef int cmp_val
        cdef int cur_val = 0
        cdef int match[6]
        cdef int cur_ctr
        cdef int dis_flag

        # reset the match array
        match[0] = match[1] = match[2] = match[3] = match[4] = match[5] = -10

        # Compare the 6-base region to AGGAGG
        limit = min(6, start - 4 - pos)
        for i in range(limit):
            if i%3 == 0 and _is_a(self.digits, self.slen, pos+i, strand):
                match[i] = 2
            elif i%3 != 0 and _is_g(self.digits, self.slen, pos+i, strand):
                match[i] = 3
            else:
                match[i] = -10

        # Find the maximally scoring motif
        max_val = 0
        for i in range(limit, 2, -1):
            for j in range(limit+1-i):
                cur_ctr = -2
                mism = 0;
                for k in range(j, j+i):
                    cur_ctr += match[k]
                    if match[k] < 0:
                        mism += 1
                if mism > 0:
                    continue
                rdis = start - (pos + j + i)
                if rdis < 5 and i < 5:
                    dis_flag = 2
                elif rdis < 5 and i >= 5:
                    dis_flag = 1
                elif rdis > 10 and rdis <= 12 and i < 5:
                    dis_flag = 1
                elif rdis > 10 and rdis <= 12 and i >= 5:
                    dis_flag = 2
                elif rdis >= 13:
                    dis_flag = 3
                else:
                    dis_flag = 0
                if rdis > 15 or cur_ctr < 6.0:
                    continue

                # Exact-Matching RBS Motifs
                if cur_ctr < 6:
                    cur_val = 0
                elif cur_ctr == 6 and dis_flag == 2:
                    cur_val = 1
                elif cur_ctr == 6 and dis_flag == 3:
                    cur_val = 2
                elif cur_ctr == 8 and dis_flag == 3:
                    cur_val = 3
                elif cur_ctr == 9 and dis_flag == 3:
                    cur_val = 3
                elif cur_ctr == 6 and dis_flag == 1:
                    cur_val = 6
                elif cur_ctr == 11 and dis_flag == 3:
                    cur_val = 10
                elif cur_ctr == 12 and dis_flag == 3:
                    cur_val = 10
                elif cur_ctr == 14 and dis_flag == 3:
                    cur_val = 10
                elif cur_ctr == 8 and dis_flag == 2:
                    cur_val = 11
                elif cur_ctr == 9 and dis_flag == 2:
                    cur_val = 11
                elif cur_ctr == 8 and dis_flag == 1:
                    cur_val = 12
                elif cur_ctr == 9 and dis_flag == 1:
                    cur_val = 12
                elif cur_ctr == 6 and dis_flag == 0:
                    cur_val = 13
                elif cur_ctr == 8 and dis_flag == 0:
                    cur_val = 15
                elif cur_ctr == 9 and dis_flag == 0:
                    cur_val = 16
                elif cur_ctr == 11 and dis_flag == 2:
                    cur_val = 20
                elif cur_ctr == 11 and dis_flag == 1:
                    cur_val = 21
                elif cur_ctr == 11 and dis_flag == 0:
                    cur_val = 22
                elif cur_ctr == 12 and dis_flag == 2:
                    cur_val = 20
                elif cur_ctr == 12 and dis_flag == 1:
                    cur_val = 23
                elif cur_ctr == 12 and dis_flag == 0:
                    cur_val = 24
                elif cur_ctr == 14 and dis_flag == 2:
                    cur_val = 25
                elif cur_ctr == 14 and dis_flag == 1:
                    cur_val = 26
                elif cur_ctr == 14 and dis_flag == 0:
                    cur_val = 27

                if tinf.rbs_wt[cur_val] < tinf.rbs_wt[max_val]:
                    continue
                if tinf.rbs_wt[cur_val] == tinf.rbs_wt[max_val] and cur_val < max_val:
                    continue
                max_val = cur_val

        return max_val

    cdef int _shine_dalgarno_mm(
        self,
        const int pos,
        const int start,
        const _training* tinf,
        const int strand
    ) nogil:
        cdef int i
        cdef int j
        cdef int k
        cdef int mism
        cdef int rdis
        cdef int limit
        cdef int max_val
        cdef int cmp_val
        cdef int cur_val = 0
        cdef int match[6]
        cdef int cur_ctr
        cdef int dis_flag

        # reset the match array
        match[0] = match[1] = match[2] = match[3] = match[4] = match[5] = -10

        # Compare the 6-base region to AGGAGG
        limit = min(6, start - 4 - pos)
        for i in range(limit):
            if i%3 == 0:
                match[i] = -3 + 5*_is_a(self.digits, self.slen, pos+i, strand)
            else:
                match[i] = -2 + 5*_is_g(self.digits, self.slen, pos+i, strand)

        # Find the maximally scoring motif
        max_val = 0
        for i in range(limit, 4, -1):
            for j in range(limit+1-i):
                cur_ctr = -2
                mism = 0;
                for k in range(j, j+i):
                    cur_ctr += match[k]
                    if match[k] < 0.0:
                        mism += 1
                        if k <= j+1 or k >= j+i-2:
                          cur_ctr -= 10
                if mism != 1:
                    continue
                rdis = start - (pos + j + i)
                if rdis < 5:
                    dis_flag = 1
                elif rdis > 10 and rdis <= 12:
                    dis_flag = 2
                elif rdis >= 13:
                    dis_flag = 3
                else:
                    dis_flag = 0
                if rdis > 15 or cur_ctr < 6:
                    continue

                # Single-Matching RBS Motifs
                if cur_ctr < 6:
                    cur_val = 0
                elif cur_ctr == 6 and dis_flag == 3:
                    cur_val = 2
                elif cur_ctr == 7 and dis_flag == 3:
                    cur_val = 2
                elif cur_ctr == 9 and dis_flag == 3:
                    cur_val = 3
                elif cur_ctr == 6 and dis_flag == 2:
                    cur_val = 4
                elif cur_ctr == 6 and dis_flag == 1:
                    cur_val = 5
                elif cur_ctr == 6 and dis_flag == 0:
                    cur_val = 9
                elif cur_ctr == 7 and dis_flag == 2:
                    cur_val = 7
                elif cur_ctr == 7 and dis_flag == 1:
                    cur_val = 8
                elif cur_ctr == 7 and dis_flag == 0:
                    cur_val = 14
                elif cur_ctr == 9 and dis_flag == 2:
                    cur_val = 17
                elif cur_ctr == 9 and dis_flag == 1:
                    cur_val = 18
                elif cur_ctr == 9 and dis_flag == 0:
                    cur_val = 19

                if tinf.rbs_wt[cur_val] < tinf.rbs_wt[max_val]:
                    continue
                if tinf.rbs_wt[cur_val] == tinf.rbs_wt[max_val] and cur_val < max_val:
                    continue
                max_val = cur_val

        return max_val

    # --- Python interface ---------------------------------------------------

    cpdef int shine_dalgarno(
        self,
        int pos,
        int start,
        TrainingInfo training_info,
        int strand=1,
        bint exact=True
    ) except -1:
        """shine_dalgarno(self, pos, start, training_info, strand=1, exact=True)\n--

        Find the highest scoring Shine-Dalgarno motif upstream of ``start``.

        Arguments:
            pos (`int`): The position where to look for the Shine-Dalgarno
                motif. Must be upstream of ``start`` (before or after,
                depending on the strand).
            start (`int`): The position of the start codon being considered.
            training_info (`~pyrodigal.TrainingInfo`): The training info
                containing the weights for the different ribosome weights.

        Keyword Arguments:
            strand (`int`): The strand to scan.
            exact (`bool`): `True` to score Shine-Dalgarno motifs matching
                exactly *AGGAGG*, `False` to allow one base mismatch.

        Returns:
            `int`: The index of the highest scoring Shine-Dalgarno motif.

        Raises:
            `ValueError`: On invalid `strand`, `pos` or `start` values.

        """
        if strand != 1 and strand != -1:
            raise ValueError(f"Invalid strand: {strand!r} (must be +1 or -1)")
        if pos < 0:
            raise ValueError(f"`pos` must be positive")
        if start < 0:
            raise ValueError(f"`start` must be positive")

        if strand == 1 and pos > start - 5:
            raise ValueError(f"`pos` is too close to `start` (must be at most `start` - 5)")
        elif strand == -1 and pos < start + 6:
            raise ValueError(f"`pos` is too close to `start` (must be at most `start + 6`)")

        cdef int phase
        with nogil:
            if exact:
                phase = self._shine_dalgarno_exact(pos, start, training_info.tinf, strand)
            else:
                phase = self._shine_dalgarno_mm(pos, start, training_info.tinf, strand)

        return phase


# --- Connection Scorer ------------------------------------------------------

_TARGET_CPU           = TARGET_CPU
_AVX2_RUNTIME_SUPPORT = False
_NEON_RUNTIME_SUPPORT = False
_SSE2_RUNTIME_SUPPORT = False
_AVX2_BUILD_SUPPORT   = False
_NEON_BUILD_SUPPORT   = False
_SSE2_BUILD_SUPPORT   = False

IF TARGET_CPU == "x86":
    cdef X86Info cpu_info = GetX86Info()
    _SSE2_RUNTIME_SUPPORT = cpu_info.features.sse2 != 0
    _AVX2_RUNTIME_SUPPORT = cpu_info.features.avx2 != 0
    _SSE2_BUILD_SUPPORT   = SSE2_BUILD_SUPPORT
    _AVX2_BUILD_SUPPORT   = AVX2_BUILD_SUPPORT
ELIF TARGET_CPU == "arm":
    cdef ArmInfo cpu_info = GetArmInfo()
    _NEON_RUNTIME_SUPPORT = cpu_info.features.neon != 0
    _NEON_BUILD_SUPPORT   = NEON_BUILD_SUPPORT
ELIF TARGET_CPU == "aarch64":
    _NEON_RUNTIME_SUPPORT = True
    _NEON_BUILD_SUPPORT   = NEON_BUILD_SUPPORT

cdef enum simd_backend:
    NONE = 0
    SSE2 = 1
    AVX2 = 2
    NEON = 3
    GENERIC = 4

cdef class ConnectionScorer:

    # --- Magic methods ------------------------------------------------------

    def __cinit__(self):
        self.capacity = 0
        self.skip_connection = self.skip_connection_raw = NULL
        self.node_types      = self.node_types_raw      = NULL
        self.node_strands    = self.node_strands_raw    = NULL
        self.node_frames     = self.node_frames_raw     = NULL

    def __init__(self, str backend="detect"):
        IF TARGET_CPU == "x86":
            if backend =="detect":
                self.backend = simd_backend.NONE
                IF SSE2_BUILD_SUPPORT:
                    if _SSE2_RUNTIME_SUPPORT:
                        self.backend = simd_backend.SSE2
                IF AVX2_BUILD_SUPPORT:
                    if _AVX2_RUNTIME_SUPPORT:
                        self.backend = simd_backend.AVX2
            elif backend == "sse":
                IF not SSE2_BUILD_SUPPORT:
                    raise RuntimeError("Extension was compiled without SSE2 support")
                ELSE:
                    if not _SSE2_RUNTIME_SUPPORT:
                        raise RuntimeError("Cannot run SSE2 instructions on this machine")
                    self.backend = simd_backend.SSE2
            elif backend == "avx":
                IF not AVX2_BUILD_SUPPORT:
                    raise RuntimeError("Extension was compiled without AVX2 support")
                ELSE:
                    if not _AVX2_RUNTIME_SUPPORT:
                        raise RuntimeError("Cannot run AVX2 instructions on this machine")
                    self.backend = simd_backend.AVX2
            elif backend == "generic":
                self.backend = simd_backend.GENERIC
            elif backend is None:
                self.backend = simd_backend.NONE
            else:
                raise ValueError(f"Unsupported backend on this architecture: {backend}")
        ELIF TARGET_CPU == "arm" or TARGET_CPU == "aarch64":
            if backend =="detect":
                self.backend = simd_backend.NONE
                IF NEON_BUILD_SUPPORT:
                    if _NEON_RUNTIME_SUPPORT:
                        self.backend = simd_backend.NEON
            elif backend == "neon":
                IF not NEON_BUILD_SUPPORT:
                    raise RuntimeError("Extension was compiled without NEON support")
                ELSE:
                    if not _NEON_RUNTIME_SUPPORT:
                        raise RuntimeError("Cannot run NEON instructions on this machine")
                    self.backend = simd_backend.NEON
            elif backend == "generic":
                self.backend = simd_backend.GENERIC
            elif backend is None:
                self.backend = simd_backend.NONE
            else:
                raise ValueError(f"Unsupported backend on this architecture: {backend}")
        ELSE:
            if backend == "detect":
                self.backend = simd_backend.NONE
            if backend == "generic":
                self.backend = simd_backend.GENERIC
            elif backend is None:
                self.backend = simd_backend.NONE
            else:
                raise ValueError(f"Unsupported backend on this architecture: {backend}")

    def __dealloc__(self):
        PyMem_Free(self.node_types_raw)
        PyMem_Free(self.node_strands_raw)
        PyMem_Free(self.node_frames_raw)
        PyMem_Free(self.skip_connection_raw)

    # --- C interface --------------------------------------------------------

    cdef int _index(self, Nodes nodes) nogil except -1:
        cdef size_t i
        # nothing to be done if we are using the Prodigal code
        if self.backend == simd_backend.NONE:
            return 0
        # reallocate if needed
        if self.capacity < nodes.length:
            with gil:
                # reallocate new memory
                self.skip_connection_raw = <uint8_t*> PyMem_Realloc(self.skip_connection_raw, nodes.length * sizeof(uint8_t) + 0x1F)
                self.node_types_raw      = <uint8_t*> PyMem_Realloc(self.node_types_raw, nodes.length      * sizeof(uint8_t) + 0x1F)
                self.node_strands_raw    = <int8_t*>  PyMem_Realloc(self.node_strands_raw, nodes.length    * sizeof(int8_t)  + 0x1F)
                self.node_frames_raw     = <uint8_t*> PyMem_Realloc(self.node_frames_raw, nodes.length     * sizeof(uint8_t) + 0x1F)
                # check that allocations were successful
                if self.skip_connection_raw == NULL:
                    raise MemoryError("Failed to allocate memory for scoring bypass index")
                if self.node_types_raw == NULL:
                    raise MemoryError("Failed to allocate memory for node type array")
                if self.node_strands_raw == NULL:
                    raise MemoryError("Failed to allocate memory for node strand array")
                if self.node_frames_raw == NULL:
                    raise MemoryError("Failed to allocate memory for node frame array")
            # record new capacity
            self.capacity = nodes.length
            # compute pointers to aligned memory
            self.skip_connection = <uint8_t*> ((<uintptr_t> self.skip_connection_raw + 0x1F) & (~0x1F))
            self.node_types      = <uint8_t*> ((<uintptr_t> self.node_types_raw      + 0x1F) & (~0x1F))
            self.node_strands    = <int8_t*>  ((<uintptr_t> self.node_strands_raw    + 0x1F) & (~0x1F))
            self.node_frames     = <uint8_t*> ((<uintptr_t> self.node_frames_raw     + 0x1F) & (~0x1F))
        # copy data from the array of nodes
        for i in range(nodes.length):
            self.node_types[i]      = nodes.nodes[i].type
            self.node_strands[i]    = nodes.nodes[i].strand
            self.node_frames[i]     = nodes.nodes[i].ndx % 3
            self.skip_connection[i] = False
        # return 0 if no exceptions were raised
        return 0

    cdef int _compute_skippable(
        self,
        int min,
        int i
    ) nogil:
        IF AVX2_BUILD_SUPPORT:
            if self.backend == simd_backend.AVX2:
                skippable_avx(self.node_strands, self.node_types, self.node_frames, min, i, self.skip_connection)
                return 0
        IF SSE2_BUILD_SUPPORT:
            if self.backend == simd_backend.SSE2:
                skippable_sse(self.node_strands, self.node_types, self.node_frames, min, i, self.skip_connection)
                return 0
        IF NEON_BUILD_SUPPORT:
            if self.backend == simd_backend.NEON:
                skippable_neon(self.node_strands, self.node_types, self.node_frames, min, i, self.skip_connection)
                return 0
        if self.backend == simd_backend.GENERIC:
            skippable_generic(self.node_strands, self.node_types, self.node_frames, min, i, self.skip_connection)
            return 0
        return 0

    @staticmethod
    cdef void _score_connection(
        Nodes nodes,
        const int p1,
        const int p2,
        const _training* tinf,
        const bint final
    ) nogil:
        cdef _node* n1 = &nodes.nodes[p1]
        cdef _node* n2 = &nodes.nodes[p2]
        cdef _node* n3

        cdef int i
        cdef int bnd
        cdef int ovlp  = 0
        cdef int maxfr = -1
        cdef int left  = n1.ndx
        cdef int right = n2.ndx

        cdef double maxval
        cdef double score   = 0.0
        cdef double scr_mod = 0.0

        # NOTE(@althonos): We can skip checking for invalid connections here
        #                  because the connection scorer already checked for
        #                  them using the SIMD filter

        # --- Edge Artifacts ---
        if n1.traceb == -1 and n1.strand == 1 and n1.type == node_type.STOP:
            return
        elif n1.traceb == -1 and n1.strand == -1 and n1.type != node_type.STOP:
            return

        # --- Genes ---
        # 5'fwd->3'fwd
        elif n1.strand == 1 and n2.strand == 1 and n1.type != node_type.STOP and n2.type == node_type.STOP:
            if n2.stop_val >= n1.ndx:
                return
            right += 2
            if final:
                score = n1.cscore + n1.sscore
            else:
                scr_mod = tinf.bias[0]*n1.gc_score[0] + tinf.bias[1]*n1.gc_score[1] + tinf.bias[2]*n1.gc_score[2]
        # 3'rev->5'rev */
        elif n1.strand == -1 and n2.strand == -1 and n1.type == node_type.STOP and n2.type != node_type.STOP:
            if n1.stop_val <= n2.ndx:
                return
            left -= 2;
            if final == 0:
                scr_mod = tinf.bias[0]*n2.gc_score[0] + tinf.bias[1]*n2.gc_score[1] + tinf.bias[2]*n2.gc_score[2]
            elif final: