plan7.pyx

# coding: utf-8
# cython: language_level=3, linetrace=True
"""High-level interface to the Plan7 data model.

Plan7 is the model architecture used by HMMER since HMMER2.

See Also:
    Details about the Plan 7 architecture in the `HMMER documentation
    <http://www.csb.yale.edu/userguides/seq/hmmer/docs/node11.html>`_.

"""

# --- C imports --------------------------------------------------------------

cimport cython
from cpython.bytes cimport PyBytes_FromStringAndSize
from cpython.list cimport PyList_New, PyList_SET_ITEM
from cpython.ref cimport PyObject
from cpython.exc cimport PyErr_Clear
from cpython.unicode cimport PyUnicode_DecodeASCII
from libc.math cimport exp, ceil
from libc.stddef cimport ptrdiff_t
from libc.stdio cimport printf
from libc.stdlib cimport calloc, malloc, realloc, free
from libc.stdint cimport uint8_t, uint32_t, int64_t
from libc.stdio cimport fprintf, FILE, stdout, fclose
from libc.string cimport memset, memcpy, memmove, strdup, strndup, strlen, strcmp, strncpy
from libc.time cimport ctime, strftime, time, time_t, tm, localtime_r
from unicode cimport PyUnicode_DATA, PyUnicode_KIND, PyUnicode_READ, PyUnicode_READY, PyUnicode_GET_LENGTH

cimport libeasel
cimport libeasel.sq
cimport libeasel.alphabet
cimport libeasel.dmatrix
cimport libeasel.fileparser
cimport libeasel.random
cimport libeasel.scorematrix
cimport libeasel.getopts
cimport libeasel.vec
cimport libhmmer
cimport libhmmer.modelconfig
cimport libhmmer.modelstats
cimport libhmmer.p7_alidisplay
cimport libhmmer.p7_hmm
cimport libhmmer.p7_builder
cimport libhmmer.p7_bg
cimport libhmmer.p7_domaindef
cimport libhmmer.p7_hit
cimport libhmmer.p7_hmmfile
cimport libhmmer.p7_pipeline
cimport libhmmer.p7_prior
cimport libhmmer.p7_profile
cimport libhmmer.p7_scoredata
cimport libhmmer.p7_tophits
cimport libhmmer.p7_trace
cimport libhmmer.tracealign
from libeasel cimport eslERRBUFSIZE, eslCONST_LOG2R
from libeasel.alphabet cimport ESL_ALPHABET, esl_alphabet_Create, esl_abc_ValidateType
from libeasel.getopts cimport ESL_GETOPTS, ESL_OPTIONS
from libeasel.sq cimport ESL_SQ
from libeasel.keyhash cimport ESL_KEYHASH
from libeasel.fileparser cimport ESL_FILEPARSER
from libhmmer cimport p7_LOCAL, p7_EVPARAM_UNSET, p7_CUTOFF_UNSET, p7_NEVPARAM, p7_NCUTOFFS, p7_offsets_e, p7_cutoffs_e, p7_evparams_e
from libhmmer.logsum cimport p7_FLogsumInit
from libhmmer.p7_builder cimport P7_BUILDER, p7_archchoice_e, p7_wgtchoice_e, p7_effnchoice_e
from libhmmer.p7_hmm cimport p7H_NTRANSITIONS, p7H_TC, p7H_GA, p7H_NC, p7H_MAP
from libhmmer.p7_hmmfile cimport p7_hmmfile_formats_e
from libhmmer.p7_hit cimport p7_hitflags_e, P7_HIT
from libhmmer.p7_alidisplay cimport P7_ALIDISPLAY
from libhmmer.p7_pipeline cimport P7_PIPELINE, p7_pipemodes_e, p7_zsetby_e, p7_strands_e, p7_complementarity_e
from libhmmer.p7_profile cimport p7_LOCAL, p7_GLOCAL, p7_UNILOCAL, p7_UNIGLOCAL
from libhmmer.p7_trace cimport P7_TRACE, p7t_statetype_e

IF HMMER_IMPL == "VMX":
    from libhmmer.impl_vmx cimport p7_oprofile, p7_omx, impl_Init
    from libhmmer.impl_vmx.io cimport p7_oprofile_Write, p7_oprofile_ReadMSV, p7_oprofile_ReadRest
    from libhmmer.impl_vmx.p7_oprofile cimport (
        P7_OPROFILE,
        p7O_NQB,
        p7_oprofile_Compare,
        p7_oprofile_Dump,
        p7_oprofile_Sizeof,
    )
ELIF HMMER_IMPL == "SSE":
    from libhmmer.impl_sse cimport p7_oprofile, p7_omx, impl_Init, p7_SSVFilter, p7O_EXTRA_SB
    from libhmmer.impl_sse.io cimport p7_oprofile_Write, p7_oprofile_ReadMSV, p7_oprofile_ReadRest
    from libhmmer.impl_sse.p7_oprofile cimport (
        P7_OPROFILE,
        p7O_NQB,
        p7_oprofile_Compare,
        p7_oprofile_Dump,
        p7_oprofile_Sizeof,
    )

from .easel cimport (
    Alphabet,
    Sequence,
    DigitalSequence,
    KeyHash,
    MSA,
    TextMSA,
    DigitalMSA,
    VectorF,
    VectorU8,
    MatrixF,
    MatrixU8,
)
from .reexports.p7_tophits cimport p7_tophits_Reuse
from .reexports.p7_hmmfile cimport (
    read_asc20hmm,
    read_asc30hmm,
    read_bin30hmm,
    v3a_magic,
    v3b_magic,
    v3c_magic,
    v3d_magic,
    v3e_magic,
    v3f_magic
)

IF UNAME_SYSNAME == "Linux":
    include "fileobj/linux.pxi"
ELIF UNAME_SYSNAME == "Darwin" or UNAME_SYSNAME.endswith("BSD"):
    include "fileobj/bsd.pxi"

# --- Python imports ---------------------------------------------------------

import array
import collections.abc
import datetime
import errno
import math
import io
import itertools
import os
import sys
import warnings

from .errors import AllocationError, UnexpectedError, AlphabetMismatch
from .utils import peekable


# --- Constants --------------------------------------------------------------

cdef dict BUILDER_ARCHITECTURE_STRATEGY = {
    "fast": p7_archchoice_e.p7_ARCH_FAST,
    "hand": p7_archchoice_e.p7_ARCH_HAND,
}

cdef dict BUILDER_WEIGHTING_STRATEGY = {
    "pb": p7_wgtchoice_e.p7_WGT_PB,
    "gsc": p7_wgtchoice_e.p7_WGT_GSC,
    "blosum": p7_wgtchoice_e.p7_WGT_BLOSUM,
    "none": p7_wgtchoice_e.p7_WGT_NONE,
    "given": p7_wgtchoice_e.p7_WGT_GIVEN,
}

cdef dict BUILDER_EFFECTIVE_STRATEGY = {
    "entropy": p7_effnchoice_e.p7_EFFN_ENTROPY,
    "exp": p7_effnchoice_e.p7_EFFN_ENTROPY_EXP,
    "clust": p7_effnchoice_e.p7_EFFN_CLUST,
    "none": p7_effnchoice_e.p7_EFFN_NONE,
}

cdef dict HMM_FILE_FORMATS = {
    "2.0": p7_hmmfile_formats_e.p7_HMMFILE_20,
    "3/a": p7_hmmfile_formats_e.p7_HMMFILE_3a,
    "3/b": p7_hmmfile_formats_e.p7_HMMFILE_3b,
    "3/c": p7_hmmfile_formats_e.p7_HMMFILE_3c,
    "3/d": p7_hmmfile_formats_e.p7_HMMFILE_3d,
    "3/e": p7_hmmfile_formats_e.p7_HMMFILE_3e,
    "3/f": p7_hmmfile_formats_e.p7_HMMFILE_3f,
}

cdef dict HMM_FILE_MAGIC = {
    v3a_magic: p7_hmmfile_formats_e.p7_HMMFILE_3a,
    v3b_magic: p7_hmmfile_formats_e.p7_HMMFILE_3b,
    v3c_magic: p7_hmmfile_formats_e.p7_HMMFILE_3c,
    v3d_magic: p7_hmmfile_formats_e.p7_HMMFILE_3d,
    v3e_magic: p7_hmmfile_formats_e.p7_HMMFILE_3e,
    v3f_magic: p7_hmmfile_formats_e.p7_HMMFILE_3f,
}

cdef dict PIPELINE_BIT_CUTOFFS = {
    "gathering": libhmmer.p7_hmm.p7H_GA,
    "noise": libhmmer.p7_hmm.p7H_NC,
    "trusted": libhmmer.p7_hmm.p7H_TC,
}

# --- Cython classes ---------------------------------------------------------


cdef class Alignment:
    """An alignment of a sequence to a profile.

    Attributes:
        domain (`Domain`): The domain this alignment corresponds to.

    """

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

    def __cinit__(self, Domain domain):
        self.domain = domain
        self._ad = domain._dom.ad

    def __len__(self):
        return self.hmm_to - self.hmm_from

    def __str__(self):
        assert self._ad != NULL

        cdef int    status
        cdef object buffer = io.BytesIO()
        cdef FILE*  fp     = fopen_obj(buffer, "w")

        try:
            status = libhmmer.p7_alidisplay.p7_nontranslated_alidisplay_Print(
                fp,
                self._ad,
                0,
                -1,
                False,
            )
            if status == libeasel.eslEWRITE:
                raise OSError("Failed to write alignment")
            elif status != libeasel.eslOK:
                raise UnexpectedError(status, "p7_alidisplay_Print")
        finally:
            fclose(fp)

        return buffer.getvalue().decode("ascii")

    cpdef VectorU8 __getstate__(self):
        cdef int       status
        cdef uint32_t  n      = 0
        cdef uint32_t  nalloc = 0
        cdef uint8_t*  buffer = NULL
        cdef VectorU8  vec

        with nogil:
            status = libhmmer.p7_alidisplay.p7_alidisplay_Serialize(self._ad, &buffer, &n, &nalloc)
        if status != libeasel.eslOK:
            raise UnexpectedError(status, "p7_alidisplay_Serialize")

        vec = VectorU8.__new__(VectorU8)
        vec._owner = None
        vec._n = vec._shape[0] = nalloc
        vec._data = <void*> buffer
        return vec

    cpdef object __setstate__(self, uint8_t[::1] state):
        cdef int      status
        cdef uint32_t offset = 0

        if self._ad == NULL:
            self.domain._dom.ad = self._ad = libhmmer.p7_alidisplay.p7_alidisplay_Create_empty()
            if self._ad == NULL:
                raise AllocationError("P7_ALIDISPLAY", sizeof(P7_ALIDISPLAY))

        with nogil:
            status = libhmmer.p7_alidisplay.p7_alidisplay_Deserialize(&state[0], &offset, self._ad)
        if status != libeasel.eslOK:
            raise UnexpectedError(status, "p7_alidisplay_Deserialize")

    # --- Properties ---------------------------------------------------------

    @property
    def hmm_from(self):
        """`int`: The start coordinate of the alignment in the query HMM.
        """
        assert self._ad != NULL
        return self._ad.hmmfrom

    @property
    def hmm_to(self):
        """`int`: The end coordinate of the alignment in the query HMM.
        """
        assert self._ad != NULL
        return self._ad.hmmto

    @property
    def hmm_name(self):
        """`bytes`: The name of the query HMM.
        """
        assert self._ad != NULL
        assert self._ad.hmmname != NULL
        return <bytes> self._ad.hmmname

    @property
    def hmm_accession(self):
        """`bytes`: The accession of the query, or its name if it has none.

        .. versionadded:: 0.1.4

        """
        assert self._ad != NULL
        assert self._ad.hmmacc != NULL
        return <bytes> self._ad.hmmacc

    @property
    def hmm_sequence(self):
        """`str`: The sequence of the query HMM in the alignment.
        """
        assert self._ad != NULL
        return self._ad.model.decode('ascii')

    @property
    def target_from(self):
        """`int`: The start coordinate of the alignment in the target sequence.
        """
        assert self._ad != NULL
        return self._ad.sqfrom

    @property
    def target_name(self):
        """`bytes`: The name of the target sequence.
        """
        assert self._ad != NULL
        assert self._ad.sqname != NULL
        return <bytes> self._ad.sqname

    @property
    def target_sequence(self):
        """`str`: The sequence of the target sequence in the alignment.
        """
        assert self._ad != NULL
        return self._ad.aseq.decode('ascii')

    @property
    def target_to(self):
        """`int`: The end coordinate of the alignment in the target sequence.
        """
        assert self._ad != NULL
        return self._ad.sqto

    @property
    def identity_sequence(self):
        """`str`: The identity sequence between the query and the target.
        """
        assert self._ad != NULL
        return self._ad.mline.decode('ascii')


cdef class Background:
    """The null background model of HMMER.

    Attributes:
        alphabet (`~pyhmmer.easel.Alphabet`): The alphabet of the
            backgound model.
        uniform (`bool`): Whether or not the null model has been created
            with uniform frequencies.
        residue_frequencies (`~pyhmmer.easel.VectorF`): The null1 background
            residue frequencies.

    .. versionchanged:: 0.4.0
       Added the `residue_frequencies` attribute.

    """

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

    def __cinit__(self):
        self._bg = NULL
        self._L  = 350
        self.alphabet = None
        self.residue_frequencies = None
        self.uniform = False

    def __init__(self, Alphabet alphabet, bint uniform=False):
        """__init__(self, alphabet, uniform=False)\n--

        Create a new background model for the given ``alphabet``.

        Arguments:
          alphabet (`pyhmmer.easel.Alphabet`): The alphabet to create the
              background model with.
          uniform (`bool`): Whether or not to create the null model with
              uniform frequencies. Defaults to `False`.

        """
        # store the alphabet so it's not deallocated
        self.alphabet = alphabet
        # store whether or not the null model has uniform frequencies
        self.uniform = uniform
        # create the background profile
        with nogil:
            if uniform:
                self._bg = libhmmer.p7_bg.p7_bg_CreateUniform(alphabet._abc)
            else:
                self._bg = libhmmer.p7_bg.p7_bg_Create(alphabet._abc)
        if self._bg == NULL:
            raise AllocationError("P7_BG", sizeof(P7_BG))
        # expose the residue frequencies as the `residue_frequencies` attribute
        self.residue_frequencies = VectorF.__new__(VectorF)
        self.residue_frequencies._data = &(self._bg.f[0])
        self.residue_frequencies._owner = self
        self.residue_frequencies._n = self.residue_frequencies._shape[0] = self.alphabet.K

    def __dealloc__(self):
        libhmmer.p7_bg.p7_bg_Destroy(self._bg)

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

    # --- Properties ---------------------------------------------------------

    @property
    def L(self):
        """`int`: The mean of the null model length distribution, in residues.
        """
        assert self._bg != NULL
        return self._L

    @L.setter
    def L(self, int L):
        assert self._bg != NULL

        cdef int    status
        cdef P7_BG* bg     = self._bg

        with nogil:
            status = libhmmer.p7_bg.p7_bg_SetLength(bg, L)
        if status != libeasel.eslOK:
            raise UnexpectedError(status, "p7_bg_SetLength")
        self._L = L

    @property
    def transition_probability(self):
        r"""`float`: The null1 transition probability (:math:`\frac{L}{L+1}`).

        .. versionadded:: 0.4.0

        """
        assert self._bg != NULL
        return self._bg.p1

    @property
    def omega(self):
        """`float`: The *prior* on *null2*/*null3*.

        .. versionadded:: 0.4.0

        """
        assert self._bg != NULL
        return self._bg.omega

    @omega.setter
    def omega(self, float omega):
        assert self._bg != NULL
        self._bg.omega = omega

    # --- Methods ------------------------------------------------------------

    cpdef Background copy(self):
        """copy(self)\n--

        Create a copy of the null model with the same parameters.

        """
        assert self._bg != NULL

        cdef Background new = Background.__new__(Background)
        with nogil:
            new._bg = libhmmer.p7_bg.p7_bg_Clone(self._bg)
        if new._bg == NULL:
            raise AllocationError("P7_BG", sizeof(P7_BG))

        new.alphabet = self.alphabet
        new.uniform = self.uniform
        new.residue_frequencies = VectorF.__new__(VectorF)
        new.residue_frequencies._data = &(new._bg.f[0])
        new.residue_frequencies._owner = new
        new.residue_frequencies._n = new.residue_frequencies._shape[0] = new.alphabet.K
        return new


cdef class Builder:
    """A factory for constructing new HMMs from raw sequences.

    Attributes:
        alphabet (`~pyhmmer.easel.Alphabet`): The alphabet the builder is
            configured to use to convert sequences to HMMs.
        randomness (`~pyhmmer.easel.Randomness`): The random number generator
            being used by the builder.
        score_matrix (`str`): The name of the substitution matrix used to
            build HMMs for single sequences.
        popen (`float`): The *gap open* probability to use when building
            HMMs from single sequences.
        pextend (`float`): The *gap extend* probability to use when building
            HMMs from single sequences.

    .. versionadded:: 0.2.0

    .. versionchanged:: 0.4.2
       Added the `~Builder.randomness` attribute.

    """

    _ARCHITECTURE_STRATEGY = dict(BUILDER_ARCHITECTURE_STRATEGY)
    _WEIGHTING_STRATEGY = dict(BUILDER_WEIGHTING_STRATEGY)
    _EFFECTIVE_STRATEGY = dict(BUILDER_EFFECTIVE_STRATEGY)

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

    def __cinit__(self):
        self._bld = NULL
        self.alphabet = None

    def __init__(
        self,
        Alphabet alphabet,
        *,
        str architecture="fast",
        str weighting="pb",
        object effective_number="entropy",
        str prior_scheme="alphabet",
        float symfrac=0.5,
        float fragthresh=0.5,
        double wid=0.62,
        double esigma=45.0,
        double eid=0.62,
        int EmL=200,
        int EmN=200,
        int EvL=200,
        int EvN=200,
        int EfL=100,
        int EfN=200,
        double Eft=0.04,
        uint32_t seed=42,
        object ere=None,
        object popen=None,
        object pextend=None,
        object score_matrix=None,
        object window_length=None,
        object window_beta=None,
    ):
        """__init__(self, alphabet, *, architecture="fast", weighting="pb", effective_number="entropy", prior_scheme="alpha", symfrac=0.5, fragthresh=0.5, wid=0.62, esigma=45.0, eid=0.62, EmL=200, EmN=200, EvL=200, EvN=200, EfL=100, EfN=200, Eft=0.04, seed=42, ere=None, popen=None, pextend=None, window_length=None, window_beta=None)\n--

        Create a new sequence builder with the given configuration.

        Arguments:
            alphabet (`~pyhmmer.easel.Alphabet`): The alphabet the builder
                expects the sequences to be in.

        Keyword Arguments:
            architecture (`str`): The algorithm to use to determine the
                model architecture, either ``"fast"`` (the default), or
                ``"hand"``.
            weighting (`str`): The algorithm to use for relative sequence
                weighting, either ``"pb"`` (the default), ``"gsc"``,
                ``"blosum"``, ``"none"``, or ``"given"``.
            effective_number (`str`, `int`, or `float`): The algorithm to
                use to determine the effective sequence number, either
                ``"entropy"`` (the default), ``"exp"``, ``"clust"``, ``"none"``.
                A number can also be given directly to set the effective
                sequence number manually.
            prior_scheme (`str`): The choice of mixture Dirichlet prior when
                parameterizing  from counts, either ``"laplace"``
                or ``"alphabet"`` (the default).
            symfrac (`float`): The residue occurrence threshold for fast
                architecture determination.
            fragthresh (`float`): A threshold such that a sequence is called
                a fragment when :math:`L \\le fragthresh \times alen`.
            wid (`double`): The percent identity threshold for BLOSUM relative
                weighting.
            esigma (`float`): The minimum total relative entropy parameter
                for effective number entropy weights.
            eid (`float`): The percent identity threshold for effective
                number clustering.
            EmL (`int`): The length of sequences generated for MSV fitting.
            EmN (`int`): The number of sequences generated for MSV fitting.
            EvL (`int`): The lenght of sequences generated for Viterbi fitting.
            EvN (`int`): The number of sequences generated for Viterbi fitting.
            EfL (`int`): The lenght of sequences generated for Forward fitting.
            EfN (`int`): The number of sequences generated for Forward fitting.
            Eft (`float`): The tail mass used for Forward fitting.
            seed (`int`): The seed to use to initialize the internal random
                number generator. If ``0`` is given, an arbitrary seed will
                be chosen based on the current time.
            ere (`double`, optional): The relative entropy target for effective
                number weighting, or `None`.
            popen (`float`, optional): The *gap open* probability to use
                when building HMMs from single sequences. The default value
                depends on the alphabet: *0.02* for proteins,
                *0.03125* for nucleotides.
            pextend (`float`, optional): The *gap extend* probability to use
                when building HMMs from single sequences. Default depends on
                the alphabet: *0.4* for proteins, *0.75* for nucleotides.
            score_matrix (`str`, optional): The name of the score matrix to
                use when building HMMs from single sequences. The only
                allowed value for nucleotide alphabets is *DNA1*. For
                proteins, *PAM30*, *PAM70*, *PAM120*, *PAM240*, *BLOSUM45*,
                *BLOSUM50*, *BLOSUM62* (the default), *BLOSUM80* or
                *BLOSUM90* can be used.
            window_length (`float`, optional): The window length for
                nucleotide sequences, essentially the max expected hit length.
                *If given, takes precedence over* ``window_beta``.
            window_beta (`float`, optional): The tail mass at which window
                length is determined for nucleotide sequences.

        """
        # extract alphabet and create raw P7_BUILDER object
        self.alphabet = alphabet
        abcty = alphabet._abc.type
        with nogil:
            self._bld = libhmmer.p7_builder.p7_builder_Create(NULL, alphabet._abc)
        if self._bld == NULL:
            raise AllocationError("P7_BG", sizeof(P7_BG))

        # create a Randomness object exposing the internal RNG
        self.randomness = Randomness.__new__(Randomness)
        self.randomness._owner = self
        self.randomness._rng = self._bld.r

        # store the seed given at initialization and reseed the RNG
        self.seed = seed

        # set numeric parameters
        self._bld.symfrac = symfrac
        self._bld.fragthresh = fragthresh
        self._bld.wid = wid
        self._bld.esigma = esigma
        self._bld.eid = eid
        self._bld.EmL = EmL
        self._bld.EmN = EmN
        self._bld.EvL = EvL
        self._bld.EvN = EvN
        self._bld.EfL = EfL
        self._bld.EfN = EfN
        self._bld.Eft = Eft

        # set the architecture strategy
        self.architecture = architecture
        _arch = BUILDER_ARCHITECTURE_STRATEGY.get(architecture)
        if _arch is not None:
            self._bld.arch_strategy = _arch
        else:
            raise ValueError(f"Invalid value for 'architecture': {architecture}")

        # set the weighting strategy
        self.weighting = weighting
        _weighting = BUILDER_WEIGHTING_STRATEGY.get(weighting)
        if _weighting is not None:
            self._bld.wgt_strategy = _weighting
        else:
            raise ValueError(f"Invalid value for 'weighting': {weighting}")

        # set the effective sequence number strategy
        self.effective_number = effective_number
        if isinstance(effective_number, (int, float)):
            self._bld.effn_strategy = p7_effnchoice_e.p7_EFFN_SET
            self._bld.eset = effective_number
        elif isinstance(effective_number, str):
            _effn = BUILDER_EFFECTIVE_STRATEGY.get(effective_number)
            if _effn is not None:
                self._bld.effn_strategy = _effn
            else:
                raise ValueError(f"Invalid value for 'effective_number': {effective_number}")
        else:
            ty = type(effective_number).__name__
            raise TypeError(f"Invalid type for 'effective_number': {ty}")

        # set the RE target if given one, otherwise the default
        # is alphabet dependent
        if ere is not None:
            self._bld.re_target = ere
        elif alphabet.is_nucleotide():
            self._bld.re_target = libhmmer.p7_ETARGET_DNA
        elif alphabet.is_amino():
            self._bld.re_target = libhmmer.p7_ETARGET_AMINO
        else:
            self._bld.re_target = libhmmer.p7_ETARGET_OTHER

        # set the prior scheme
        self.prior_scheme = prior_scheme
        if prior_scheme is None:
            self._bld.prior = NULL
        elif prior_scheme == "laplace":
            self._bld.prior = libhmmer.p7_prior.p7_prior_CreateLaplace(self.alphabet._abc)
        elif prior_scheme == "alphabet":
            if alphabet.is_amino():
                self._bld.prior = libhmmer.p7_prior.p7_prior_CreateAmino()
            elif alphabet.is_nucleotide():
                self._bld.prior = libhmmer.p7_prior.p7_prior_CreateNucleic()
            else:
                self._bld.prior = libhmmer.p7_prior.p7_prior_CreateLaplace(self.alphabet._abc)
        else:
            raise ValueError("Invalid value for 'prior_scheme': {prior_scheme}")

        # set the gap probabilities and score matrix using alphabet-specific
        # defaults (this is only used when building a HMM for a single sequence)
        if alphabet.is_nucleotide():
            self.score_matrix = "DNA1" if score_matrix is None else score_matrix
            self.popen = 0.03125 if popen is None else popen
            self.pextend = 0.75 if pextend is None else pextend
        else:
            self.score_matrix = "BLOSUM62" if score_matrix is None else score_matrix
            self.popen = 0.02 if popen is None else popen
            self.pextend = 0.4 if pextend is None else pextend

        # set the window options
        self.window_length = window_length
        if window_beta is None:
            self.window_beta = libhmmer.p7_builder.p7_DEFAULT_WINDOW_BETA
        else:
            self.window_beta = window_beta

    def __dealloc__(self):
        libhmmer.p7_builder.p7_builder_Destroy(self._bld)

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

    # --- Properties ---------------------------------------------------------

    @property
    def seed(self):
        """`int`: The seed given at builder initialization.

        Setting this attribute to a different value will cause the internal
        random number generator to be reseeded immediately.

        .. versionchanged:: 0.4.2
           Avoid shadowing initial null seeds given on builder initialization.

        """
        return self._seed

    @seed.setter
    def seed(self, uint32_t seed):
        self._seed = seed
        self._bld.do_reseeding = seed != 0
        self.randomness._seed(seed)

    @property
    def window_length(self):
        """`int` or `None`: The window length for nucleotide sequences.
        """
        assert self._bld != NULL
        return None if self._bld.w_len == -1 else self._bld.w_len

    @window_length.setter
    def window_length(self, object window_length):
        assert self._bld != NULL
        if window_length is None:
            self._bld.w_len = -1
        elif window_length > 3:
            self._bld.w_len = window_length
        else:
            raise ValueError(f"Invalid window length: {window_length!r}")

    @property
    def window_beta(self):
        """`float`: The tail mass at which window length is determined.
        """
        assert self._bld != NULL
        return self._bld.w_beta

    @window_beta.setter
    def window_beta(self, double window_beta):
        assert self._bld != NULL
        if window_beta > 1 or window_beta < 0:
            raise ValueError(f"Invalid window tail mass: {window_beta!r}")
        self._bld.w_beta = window_beta

    # --- Methods ------------------------------------------------------------

    cpdef tuple build(
        self,
        DigitalSequence sequence,
        Background background,
    ):
        """build(self, sequence, background)\n--

        Build a new HMM from ``sequence`` using the builder configuration.

        Arguments:
            sequence (`~pyhmmer.easel.DigitalSequence`): A single biological
                sequence in digital mode to build a HMM with.
            background (`pyhmmer.plan7.background`): The background model
                to use to create the HMM.

        Returns:
            (`HMM`, `Profile`, `OptimizedProfile`): A tuple containing the
            new HMM as well as profiles to be used directly in a `Pipeline`.

        Raises:
            `~pyhmmer.errors.AlphabetMismatch`: When either ``sequence`` or
                ``background`` have the wrong alphabet for this builder.
            `ValueError`: When the HMM cannot be created successfully
                from the input; the error message contains more details.

        Hint:
            The score matrix and the gap probabilities used here can be set
            when initializing the builder, or changed by setting a new value
            to the right attribute::

                >>> alphabet = easel.Alphabet.amino()
                >>> background = plan7.Background(alphabet)
                >>> builder = plan7.Builder(alphabet, popen=0.05)
                >>> builder.score_matrix = "BLOSUM90"
                >>> hmm, _, _ = builder.build(proteins[0], background)

        .. versionchanged:: 0.4.6
           Sets the `HMM.creation_time` attribute with the current time.

        .. versionchanged:: 0.4.10
           The score system is now cached between calls to `Builder.build`.

        """
        assert self._bld != NULL

        cdef int              status
        cdef HMM              hmm     = HMM.__new__(HMM)
        cdef Profile          profile = Profile.__new__(Profile)
        cdef OptimizedProfile opti    = OptimizedProfile.__new__(OptimizedProfile)
        cdef bytes            mx      = self.score_matrix.encode('utf-8')
        cdef char*            mx_ptr  = <char*> mx
        cdef str              msg

        # reseed RNG used by the builder if needed
        if self._bld.do_reseeding:
            self.randomness._seed(self.seed)

        # check alphabet and assign it to newly created objects
        hmm.alphabet = profile.alphabet = opti.alphabet = self.alphabet
        if not self.alphabet._eq(background.alphabet):
            raise AlphabetMismatch(self.alphabet, background.alphabet)
        if not self.alphabet._eq(sequence.alphabet):
            raise AlphabetMismatch(self.alphabet, sequence.alphabet)

        # only load the score system if it hasn't been loaded already,
        # or if a different score system is in use
        if self._bld.S == NULL or strcmp(self._bld.S.name, mx_ptr) != 0:
            with nogil:
                status = libhmmer.p7_builder.p7_builder_LoadScoreSystem(
                    self._bld,
                    mx_ptr,
                    self.popen,
                    self.pextend,
                    background._bg
                )
            if status == libeasel.eslEINVAL:
                msg = self._bld.errbuf.decode('utf-8', 'ignore')
                raise RuntimeError(f"failed to set single sequence score system: {msg}")
            elif status != libeasel.eslOK:
                raise UnexpectedError(status, "p7_builder_LoadScoreSystem")
        else:
            self._bld.popen = self.popen
            self._bld.pextend = self.pextend

        # build HMM and profiles
        with nogil:
            status = libhmmer.p7_builder.p7_SingleBuilder(
                self._bld,
                sequence._sq,
                background._bg,
                &hmm._hmm, # HMM
                NULL, # traceback
                &profile._gm, # profile,
                &opti._om, # optimized profile
            )
        if status == libeasel.eslOK:
            hmm.command_line = " ".join(sys.argv)
        elif status == libeasel.eslEINVAL:
            msg = self._bld.errbuf.decode("utf-8", "replace")
            raise ValueError("Could not build HMM: {}".format(msg))
        else:
            raise UnexpectedError(status, "p7_SingleBuilder")

        # return newly built HMM, profile and optimized profile
        return hmm, profile, opti

    cpdef tuple build_msa(
        self,
        DigitalMSA msa,
        Background background,
    ):
        """build_msa(self, msa, background)\n--

        Build a new HMM from ``msa`` using the builder configuration.

        Arguments:
            msa (`~pyhmmer.easel.DigitalMSA`): A multiple sequence
                alignment in digital mode to build a HMM with.
            background (`pyhmmer.plan7.background`): The background model
                to use to create the HMM.

        Returns:
            (`HMM`, `Profile`, `OptimizedProfile`): A tuple containing the
            new HMM as well as profiles to be used directly in a `Pipeline`.

        Raises:
            `~pyhmmer.errors.AlphabetMismatch`: When either ``msa`` or
                ``background`` have the wrong alphabet for this builder.
            `ValueError`: When the HMM cannot be created successfully
                from the input; the error message contains more details.

        Caution:
            HMMER requires that every HMM has a name, so the `Builder` will
            attempt to use the name of the sequence to name the HMM. Passing
            an MSA without a name will result in an error::

                >>> alphabet = easel.Alphabet.amino()
                >>> msa = easel.TextMSA(sequences=[
                ...   easel.TextSequence(name=b"ustiA", sequence="YAIG"),
                ...   easel.TextSequence(name=b"ustiB", sequence="YVIG")
                ... ])
                >>> builder = plan7.Builder(alphabet)
                >>> background = plan7.Background(alphabet)
                >>> builder.build_msa(msa.digitize(alphabet), background)
                Traceback (most recent call last):
                  ...
                ValueError: Could not build HMM: Unable to name the HMM.

        .. versionadded:: 0.3.0

        .. versionchanged:: 0.4.6
           Sets the `HMM.creation_time` attribute with the current time.

        """
        assert self._bld != NULL

        cdef int              status
        cdef HMM              hmm     = HMM.__new__(HMM)
        cdef Profile          profile = Profile.__new__(Profile)
        cdef OptimizedProfile opti    = OptimizedProfile.__new__(OptimizedProfile)
        cdef str              msg
        cdef size_t           k

        # unset builder probabilities and score system if they were set
        # by a call to `build`, since we won't be needing them here
        self._bld.popen = self._bld.pextend = -1
        if self._bld.S != NULL:
            libeasel.scorematrix.esl_scorematrix_Destroy(self._bld.S)
            self._bld.S = NULL
        if self._bld.Q != NULL:
            libeasel.dmatrix.esl_dmatrix_Destroy(self._bld.Q)
            self._bld.Q = NULL

        # reseed RNG used by the builder if needed
        if self._bld.do_reseeding:
            self.randomness._seed(self.seed)

        # check alphabet and assign it to newly created objects
        hmm.alphabet = profile.alphabet = opti.alphabet = self.alphabet
        if not self.alphabet._eq(background.alphabet):
            raise AlphabetMismatch(self.alphabet, background.alphabet)
        if not self.alphabet._eq(msa.alphabet):
            raise AlphabetMismatch(self.alphabet, msa.alphabet)

        # build HMM
        with nogil:
            # build HMM and profiles
            status = libhmmer.p7_builder.p7_Builder(
                self._bld,
                msa._msa,
                background._bg,
                &hmm._hmm, # HMM
                NULL, # traceback
                &profile._gm, # profile,
                &opti._om, # optimized profile
                NULL # modified msa
            )
        if status == libeasel.eslOK:
            hmm.command_line = " ".join(sys.argv)
        elif status == libeasel.eslEINVAL:
            msg = self._bld.errbuf.decode("utf-8", "replace")
            raise ValueError("Could not build HMM: {}".format(msg))
        else:
            raise UnexpectedError(status, "p7_Builder")

        # return newly built HMM, profile and optimized profile
        return hmm, profile, opti

    cpdef Builder copy(self):
        """copy(self)\n--

        Create a duplicate `Builder` instance with the same arguments.

        """
        assert self._bld != NULL
        return Builder(
            self.alphabet,
            architecture=self.architecture,
            weighting=self.weighting,
            effective_number=self.effective_number,
            prior_scheme=self.prior_scheme,
            symfrac=self._bld.symfrac,
            fragthresh=self._bld.fragthresh,
            wid=self._bld.wid,
            esigma=self._bld.esigma,
            eid=self._bld.eid,
            EmL=self._bld.EmL,
            EmN=self._bld.EmN,
            EvL=self._bld.EvL,
            EvN=self._bld.EvN,
            EfL=self._bld.EfL,
            EfN=self._bld.EfN,