aspartik.b3

class Prior:

@runtime_checkable
class Prior(Protocol):
    """
    Interface which describes all prior distributions

    A `Prior` object will be queried by `MCMC` on each step after the operator
    edits the state to get the prior probability of the new state.  It's the
    responsibility of the prior to track the stateful objects it's interested
    in, so they will typically be passed in the constructor.

    Variable priors, coalescents, birth-death models, and all other
    non-likelihood models implement this protocol.
    """

    def probability(self) -> float:
        """Calculates the log prior probability of the model state

        The return value must be a **natural logarithm** of the probability.

        `MCMC` will short-circuit and abort the move if `probability` returns a
        negative infinity.  This can be used to avoid expensive likelihood
        calculations for obviously invalid moves, like going out of variable
        bounds.
        """
        ...
#

Interface which describes all prior distributions

A Prior object will be queried by MCMC on each step after the operator edits the state to get the prior probability of the new state. It's the responsibility of the prior to track the stateful objects it's interested in, so they will typically be passed in the constructor.

Variable priors, coalescents, birth-death models, and all other non-likelihood models implement this protocol.

def probability(self) -> float

    def probability(self) -> float:
        """Calculates the log prior probability of the model state

        The return value must be a **natural logarithm** of the probability.

        `MCMC` will short-circuit and abort the move if `probability` returns a
        negative infinity.  This can be used to avoid expensive likelihood
        calculations for obviously invalid moves, like going out of variable
        bounds.
        """
        ...
#

Calculates the log prior probability of the model state

The return value must be a natural logarithm of the probability.

MCMC will short-circuit and abort the move if probability returns a negative infinity. This can be used to avoid expensive likelihood calculations for obviously invalid moves, like going out of variable bounds.

class Operator:

class Operator(Protocol):
    """
    Objects which propose moves by editing state

    It is the responsibility of the objects to track the parts of the state it
    might want to edit.  Typically these objects will be passed in the
    constructor.
    """

    weight: float
    """Influences the probability of the operator being picked

    On each step `MCMC` picks a random operator from the list passed to it.
    It uses this value to weight them.  So, the larger it is, the more
    often the operator will be picked, and visa versa.  This value is read
    once on startup.  Therefore, if it's changed mid-execution the old
    cached value will still be used.
    """

    def propose(self) -> Proposal:
        """Proposes a new MCMC step

        It is presumed that the operator will store all the references to
        parameters and trees it wants to edit and will change them accordingly.
        If a move cannot be proposed for any reason `Proposal.Reject` should be
        returned.  `MCMC` will deal with rolling back the state.
        """
        ...
#

Objects which propose moves by editing state

It is the responsibility of the objects to track the parts of the state it might want to edit. Typically these objects will be passed in the constructor.

weight: float

#

Influences the probability of the operator being picked

On each step MCMC picks a random operator from the list passed to it. It uses this value to weight them. So, the larger it is, the more often the operator will be picked, and visa versa. This value is read once on startup. Therefore, if it's changed mid-execution the old cached value will still be used.

def propose(self) -> aspartik.b3.Proposal

    def propose(self) -> Proposal:
        """Proposes a new MCMC step

        It is presumed that the operator will store all the references to
        parameters and trees it wants to edit and will change them accordingly.
        If a move cannot be proposed for any reason `Proposal.Reject` should be
        returned.  `MCMC` will deal with rolling back the state.
        """
        ...
#

Proposes a new MCMC step

It is presumed that the operator will store all the references to parameters and trees it wants to edit and will change them accordingly. If a move cannot be proposed for any reason Proposal.Reject should be returned. MCMC will deal with rolling back the state.

class Callback:

class Callback(Protocol):
    """
    Custom callbacks

    `b3` supports arbitrary logging and checks via this protocol.  The `call`
    function is passed a reference to the main `MCMC` object, so the callback
    can either take state variables in the constructor of fetch them via the
    `MCMC` attributes.

    The `call` function won't be called after each step for efficiency.  See
    [`every`](#Callback.every) for configuring how often the callback will be
    invoked.
    """

    every: int
    """How often this callback should be called

    The `MCMC` will call each callback object when `index % every` is 0.  This
    value is read once when MCMC is created, so if it's changed during
    execution, the old `every` value will continue to be used.
    """

    def call(self, mcmc: MCMC) -> None:
        """
        A custom operation

        Used by loggers and other periodic actions.
        """
#

Custom callbacks

b3 supports arbitrary logging and checks via this protocol. The call function is passed a reference to the main MCMC object, so the callback can either take state variables in the constructor of fetch them via the MCMC attributes.

The call function won't be called after each step for efficiency. See every for configuring how often the callback will be invoked.

every: int

#

How often this callback should be called

The MCMC will call each callback object when index % every is 0. This value is read once when MCMC is created, so if it's changed during execution, the old every value will continue to be used.

def call(self, mcmc: aspartik.b3.MCMC) -> None

    def call(self, mcmc: MCMC) -> None:
        """
        A custom operation

        Used by loggers and other periodic actions.
        """
#

A custom operation

Used by loggers and other periodic actions.

class Stateful:

@runtime_checkable
class Stateful(Protocol):
    """
    Epoch-versioned objects for use with MCMC

    All stateful objects handled by `MCMC` must conform to this protocol.
    During each step operators can edit objects using whatever APIs provided.
    Then, at the end of the step, `MCMC` calls `accept` if the move has been
    accepted, or `reject` otherwise.
    """

    def accept(self) -> None:
        """Accept changes made during the current step"""

    def reject(self) -> None:
        """Reject changes made during the current step

        This method must roll the state of the object back to how it was at the
        beginning of the MCMC step.
        """
#

Epoch-versioned objects for use with MCMC

All stateful objects handled by MCMC must conform to this protocol. During each step operators can edit objects using whatever APIs provided. Then, at the end of the step, MCMC calls accept if the move has been accepted, or reject otherwise.

def accept(self) -> None

    def accept(self) -> None:
        """Accept changes made during the current step"""
#

Accept changes made during the current step

def reject(self) -> None

    def reject(self) -> None:
        """Reject changes made during the current step

        This method must roll the state of the object back to how it was at the
        beginning of the MCMC step.
        """
#

Reject changes made during the current step

This method must roll the state of the object back to how it was at the beginning of the MCMC step.

class MCMC:

#

The main object which runs the analysis

callbacks

#

A list of callbacks

operators

#

A list of active operators

posterior

#

Posterior probability for the last accepted step

current_step

#

Index of the current MCMC step

Starts from 0, includes burn-in.

prior

#

Prior likelihood for the current step

Note that unlike posterior and Likelihood, this property isn't cached. It will trigger a recalculation on all priors on each access.

rng

#

Randomness source of this analysis

This objects is used for internal randomness generation (such as picking the operator on each step). Since the underlying object is shared, using it will outside of MCMC might alter the rest of the analysis.

operator_statistics

#

Operator statistics for this run

Returns a list of (operator, results, propose, likelihood) tuples for each operator. propose and likelihood records the total time the MCMC spent waiting for the operator to generate a proposal and calculate it respectively. operator is the reference to the original operator object. And results is a list of step results.

priors

#

All priors

state

#

All stateful objects tracked by this MCMC instance

likelihood

#

Likelihood calculator object

def run(self, /, n)

#

Execute n steps of the Markov chain

This yields flow control to the Rust core until the simulation is done. Press Ctrl+C to interrupt and stop the execution.

def measure_operator(self, /, operator_index, length)

#

TODO: refine and document

class Proposal:

#

A result of the move proposed by an operator

While the operators edit the tree directly, they need to communicate the status of their move to MCMC. This is the class used for that.

def Reject(cls, /)

#

Aborts the move unconditionally

All of the trees and parameters are rolled back. This is relatively fast, as it typically skips recalculating the likelihoods.

def Hastings(cls, /, ratio)

#

Proposes the move with the ratio

This is the ratio from the Metropolis–Hastings algorithm.

def Accept(cls, /)

#

Accepts the move unconditionally

class Tree:

#

A phylogenetic tree

Unlike BEAST2, where the tree is implemented as a collection of nodes pointing to each other, in b3 Tree is a self-contained data structure which holds all of the topology and heights. This means that nodes (Internal and Leaf) are identifiers of nodes in a given Tree object, much like indices of an array. This means that all operations, such as getting parents of node heights, have to go through Tree's methods.

The current implementation only supports bifurcating trees.

names

#

A list of all leaf names.

The order is the same as leaf indices: the first name is that of Leaf(0), the second one is Leaf(1), and so on.

num_leaves

#

Number of leaf nodes

num_edges

#

Total number of edges

num_nodes

#

Total number of nodes in the tree

root

#

Returns the root node of the tree

Note that the root node might change after tree has been edited, so the returned node is only guaranteed to be root as long as the tree hasn't been edited.

num_internals

#

Number of internal nodes (those with children)

def from_newick(cls, /, newick)

#

Initializes the tree from a Newick object.

The Newick tree must be strictly bifurcating and all of its edges must have a defined length.

def set_random_edges(self, /, rng)

#

Randomizes the tree structure

This methods creates a random Prüfer sequence and rearranges the graph according to it. Note that it will always the internal node with the largest index (num_nodes - 1) the root.

def set_random_heights(self, /, diff, rng)

#

Randomizes the heights of internal nodes

Each internal node gets a height distributed uniformly between diff and 2 * diff plus the height of the highest of its children.

def scale(self, /, scale)

#

Multiplies the heights of all internal nodes by scale

Exceptions

Throws a RuntimeError if any of the internal nodes would be moved below either of its children.

def update_edge(self, /, edge, new_child)

#

Sets the child of edge to new_child

This will only change the child, so the parent (internal node from which edge comes out) will now have node as a child.

This function doesn't do any validation, it's up to the operator to preserve the validity of the tree.

def set_height(self, /, node, height)

#

Sets the height of node to height

def set_root(self, /, node)

#

Makes node the root of the tree

As the topology can be temporarily broken while the edges are being swapped, Tree can't automatically figure out which node is the root one. So, operators which change the root of the tree have to update it manually.

def swap_parents(self, /, a, b)

#

Swaps the parents of nodes a and b

a and b must not be a descendant/ancestors and neither of them can be a root node. If a and b share the same parent, they switch polarity (left child becomes the right child and visa versa).

def is_internal(self, /, node)

#

Returns True if the node is internal

def is_leaf(self, /, node)

#

Returns True if the node is a leaf

def height_of(self, /, node)

#

Returns the height of node

Height here means node's age in some unlabeled units.

def children_of(self, /, node)

#

Returns a tuple of the left and right children of node

This function takes the Internal type as its input, so it is guaranteed to always return the children. See as_internal for converting general nodes to internal ones.

def other_child(self, /, parent, child)

#

Returns the child of parent other than child

Throws an error if child isn't a child of parent.

def random_intersecting_edge(self, /, height, rng)

#

Returns a random edge which intersects height

"Intersects" here means that the edge parent is higher than height and the child is lower. The comparisons are strict: if either node is exactly at height, the edge won't be picked.

Returns None if there is no such node.

def edge_index(self, /, child)

#

Returns the index of an edge from child to its parent

def edge_length(self, /, edge)

#

Returns the length of edge

The length is the distance between the parent and the child nodes of that edge

def edge_nodes(self, /, edge)

#

Returns the (child, parent) tuple corresponding to an edge

def parent_of(self, /, node)

#

Returns the parent of node, or None for the root node

def is_grandparent(self, /, node)

#

Returns True if both children of this node are also internal

def num_grandparents(self, /)

#

Number of nodes for whom is_grandparent returns True

def nodes(self, /)

#

An iterator over all of trees nodes

All of the Leaf nodes go before Internal ones.

def internals(self, /)

#

Returns an iterator over internal nodes.

def leaves(self, /)

#

Returns an iterator over all leaves.

def random_node(self, /, rng)

#

Samples a random node from the tree.

def random_nonroot_node(self, /, rng)

#

Samples a random non-root node

def random_internal(self, /, rng)

#

Samples a random internal node

def random_nonroot_internal(self, /, rng)

#

Samples a random non-root internal node

def random_leaf(self, /, rng)

#

Samples a random leaf node from a tree.

def leaf_by_name(self, /, name)

#

Gets a named leaf or None if the name is not found

def total_length(self, /)

#

The total length of all tree edges

def has_dated_tips(self, /)

#

Returns true if any of the leaves have a non-0 height

def validate(self, /)

#

Throws an exception if a tree is malformed

This function ensures that:

  • No leaf has become anyone's parent.
  • All parent nodes are older than their children.
  • Parents match their children (mismatches can happen when update_edge is used incorrectly).
  • There's only one root (two or more can be set with set_root).
  • The tree is a tree, meaning that topologically it has no cycles and is connected.

def newick(self, /, internal_ids=False)

#

Returns the tree topology in the Newick format

Leaf nodes will be labeled with the names passed to the constructor while the internal nodes are unlabeled.

def accept(self, /)

#

def reject(self, /)

#

Node

#

Any node of the phylogenetic tree

Used for type hints in places where there isn't a need to distinguish between internal and leaf nodes.