MCMC Moves

Reversible tree moves for Metropolis–Hastings sampling in tree space.

Each move returns (new_tree, log_hastings_ratio).

nni_move

nni_move(tree: HifukuTree, taxon_table: TaxonTable, rng: Generator) -> tuple[HifukuTree, float]

Nearest-neighbor interchange (NNI).

For a randomly chosen internal non-root node u, swaps one of u's children with u's sibling across the edge connecting u to its parent.

For a binary tree with n leaves there are always n-2 valid NNI edges and 2 variants each, giving 2(n-2) moves regardless of topology. The proposal is therefore symmetric and log_hastings_ratio = 0.

Parameters:

Name	Type	Description	Default
tree	HifukuTree		required
taxon_table	TaxonTable		required
rng	Generator		required

Returns:

Name	Type	Description
new_tree	HifukuTree
log_hastings_ratio	float	Always 0.0; -inf when n < 3 (no valid NNI edges).

Notes

• Minh et al. (2020) : doNNI, NNIMove struct.

Source code in src/hifuku/moves.py

def nni_move(
    tree: HifukuTree,
    taxon_table: TaxonTable,
    rng: np.random.Generator,
) -> tuple[HifukuTree, float]:
    """
    Nearest-neighbor interchange (NNI).

    For a randomly chosen internal non-root node ``u``, swaps one of
    ``u``'s children with ``u``'s sibling across the edge connecting ``u``
    to its parent.

    For a binary tree with n leaves there are always n-2 valid NNI edges
    and 2 variants each, giving 2(n-2) moves regardless of topology. The
    proposal is therefore symmetric and ``log_hastings_ratio = 0``.

    Parameters
    ----------
    tree : HifukuTree
    taxon_table : TaxonTable
    rng : np.random.Generator

    Returns
    -------
    new_tree : HifukuTree
    log_hastings_ratio : float
        Always 0.0; ``-inf`` when n < 3 (no valid NNI edges).

    Notes
    -----
    - Minh et al. (2020) : doNNI, NNIMove struct.
    """
    # All internal non-root nodes are valid: parent index >= n_leaves for all.
    # (Root index = n_nodes-1 >= n_leaves always.)
    candidates = list(range(tree.n_leaves, tree.n_nodes - 1))
    if not candidates:
        return tree, -math.inf

    u = int(rng.choice(candidates))
    p = int(tree.parent[u])

    # Sibling of u under p
    lc_p = int(tree.left_child[p])
    s = int(tree.right_child[p]) if lc_p == u else lc_p

    # Pick which child of u to swap with s
    a = int(tree.left_child[u])
    b = int(tree.right_child[u])
    swap_child = a if bool(rng.integers(0, 2)) else b

    new_parent = tree.parent.copy()
    new_lc = tree.left_child.copy()
    new_rc = tree.right_child.copy()

    # Reparent: swap_child → p, s → u
    new_parent[swap_child] = p
    new_parent[s] = u

    # Update p's child pointer (s → swap_child)
    if int(tree.left_child[p]) == s:
        new_lc[p] = swap_child
    else:
        new_rc[p] = swap_child

    # Update u's child pointer (swap_child → s)
    if int(tree.left_child[u]) == swap_child:
        new_lc[u] = s
    else:
        new_rc[u] = s

    tmp = HifukuTree(
        parent=new_parent, left_child=new_lc, right_child=new_rc,
        edge_len=tree.edge_len.copy(), names=list(tree.names),
        n_nodes=tree.n_nodes, n_leaves=tree.n_leaves, root=tree.root,
    )
    tmp.node_order = tmp.compute_node_order()
    if tree.global_leaf_indices is not None:
        tmp.global_leaf_indices = tree.global_leaf_indices.copy()
    return tmp, 0.0

spr_move

spr_move(tree: HifukuTree, taxon_table: TaxonTable, rng: Generator) -> tuple[HifukuTree, float]

Subtree prune-regraft (SPR).

Detaches a random subtree and reattaches it at a random edge in the remaining tree, splitting that edge evenly.

The proposal: (1) pick a prune edge uniformly from n_valid_prune(T) valid edges (pruned subtree has ≤ n-2 leaves); (2) pick a graft edge uniformly from the remaining tree. The graft count n_graft = 2*(n - leaf_count(prune)) - 2 cancels in the Hastings ratio:

.. code-block:: text

log_hastings = log(n_valid_prune(T)) - log(n_valid_prune(T'))

Parameters:

Name	Type	Description	Default
tree	HifukuTree		required
taxon_table	TaxonTable		required
rng	Generator		required

Returns:

Name	Type	Description
new_tree	HifukuTree
log_hastings_ratio	float	Exact Hastings ratio as above; -inf when n < 4.

Notes

• Minh et al. (2020) : SPR framework.

Source code in src/hifuku/moves.py

def spr_move(
    tree: HifukuTree,
    taxon_table: TaxonTable,
    rng: np.random.Generator,
) -> tuple[HifukuTree, float]:
    """
    Subtree prune-regraft (SPR).

    Detaches a random subtree and reattaches it at a random edge in the
    remaining tree, splitting that edge evenly.

    The proposal: (1) pick a prune edge uniformly from ``n_valid_prune(T)``
    valid edges (pruned subtree has ≤ n-2 leaves); (2) pick a graft edge
    uniformly from the remaining tree. The graft count
    ``n_graft = 2*(n - leaf_count(prune)) - 2`` cancels in the Hastings ratio:

    .. code-block:: text

        log_hastings = log(n_valid_prune(T)) - log(n_valid_prune(T'))

    Parameters
    ----------
    tree : HifukuTree
    taxon_table : TaxonTable
    rng : np.random.Generator

    Returns
    -------
    new_tree : HifukuTree
    log_hastings_ratio : float
        Exact Hastings ratio as above; ``-inf`` when n < 4.

    Notes
    -----
    - Minh et al. (2020) : SPR framework.
    """
    total_leaves = tree.n_leaves
    root_n = _from_tree(tree)

    all_e = _all_edges(root_n)
    valid_prune = [(p, c) for p, c in all_e if c.leaf_count() <= total_leaves - 2]
    if not valid_prune:
        return tree, -math.inf

    n_fwd_prune = len(valid_prune)

    # Deep copy for mutation
    root_copy = deepcopy(root_n)
    all_e_copy = _all_edges(root_copy)
    valid_copy = [(p, c) for p, c in all_e_copy if c.leaf_count() <= total_leaves - 2]

    pi = int(rng.integers(0, n_fwd_prune))
    prune_par, prune_node = valid_copy[pi]

    # Detach prune_node from its parent
    prune_par.children = [c for c in prune_par.children if c.nid != prune_node.nid]

    # Collapse prune_par if degree-1 after detachment
    if len(prune_par.children) == 1:
        survivor = prune_par.children[0]
        if prune_par.nid == root_copy.nid:
            root_copy = survivor
            root_copy.edge_len = float("nan")
        else:
            gpar = _find_parent(root_copy, prune_par.nid)
            if gpar is not None:
                survivor.edge_len = prune_par.edge_len + survivor.edge_len
                gpar.children = [
                    survivor if c.nid == prune_par.nid else c
                    for c in gpar.children
                ]

    # Collect graft edges in remaining tree and pick one
    graft_edges = _all_edges(root_copy)
    if not graft_edges:
        return tree, -math.inf

    gi = int(rng.integers(0, len(graft_edges)))
    graft_par, graft_child = graft_edges[gi]

    # Insert prune_node by splitting graft edge
    old_len = graft_child.edge_len
    half = old_len / 2.0 if math.isfinite(old_len) and old_len > 0.0 else 0.0
    new_node = _N(nid=-1, name="", edge_len=half,
                  children=[graft_child, prune_node])
    graft_child.edge_len = half
    graft_par.children = [
        new_node if c.nid == graft_child.nid else c
        for c in graft_par.children
    ]

    # Hastings ratio: only n_valid_prune differs between forward and reverse
    n_bwd_prune = _n_valid_prune(root_copy, total_leaves)
    log_hastings = math.log(n_fwd_prune) - math.log(max(1, n_bwd_prune))

    return _n_to_tree(root_copy, taxon_table), log_hastings

branch_scale_move

branch_scale_move(tree: HifukuTree, taxon_table: TaxonTable, rng: Generator, width: float = 0.5) -> tuple[HifukuTree, float]

Global branch-length scaling.

Multiplies all finite edge lengths by exp(U), U ~ Uniform(-width/2, width/2). The Jacobian of the transformation t → t·exp(U) for each of n_branches independent branches gives log_hastings_ratio = n_branches · U.

Parameters:

Name	Type	Description	Default
tree	HifukuTree		required
taxon_table	TaxonTable		required
rng	Generator		required
width	float	Half-width of the uniform proposal window (default 0.5).	0.5

Returns:

Name	Type	Description
new_tree	HifukuTree
log_hastings_ratio	float	n_branches * U (Jacobian term).

Source code in src/hifuku/moves.py

def branch_scale_move(
    tree: HifukuTree,
    taxon_table: TaxonTable,
    rng: np.random.Generator,
    width: float = 0.5,
) -> tuple[HifukuTree, float]:
    """
    Global branch-length scaling.

    Multiplies all finite edge lengths by ``exp(U)``,
    ``U ~ Uniform(-width/2, width/2)``. The Jacobian of the transformation
    ``t → t·exp(U)`` for each of ``n_branches`` independent branches gives
    ``log_hastings_ratio = n_branches · U``.

    Parameters
    ----------
    tree : HifukuTree
    taxon_table : TaxonTable
    rng : np.random.Generator
    width : float
        Half-width of the uniform proposal window (default 0.5).

    Returns
    -------
    new_tree : HifukuTree
    log_hastings_ratio : float
        ``n_branches * U`` (Jacobian term).
    """
    finite_mask = np.isfinite(tree.edge_len)
    n_branches = int(finite_mask.sum())
    if n_branches == 0:
        return tree, 0.0

    U = float(rng.uniform(-width / 2.0, width / 2.0))
    factor = math.exp(U)

    new_el = np.where(finite_mask, tree.edge_len * factor, tree.edge_len)

    tmp = HifukuTree(
        parent=tree.parent.copy(), left_child=tree.left_child.copy(),
        right_child=tree.right_child.copy(), edge_len=new_el,
        names=list(tree.names), n_nodes=tree.n_nodes,
        n_leaves=tree.n_leaves, root=tree.root,
        node_order=tree.node_order.copy() if tree.node_order is not None else None,
    )
    if tree.global_leaf_indices is not None:
        tmp.global_leaf_indices = tree.global_leaf_indices.copy()
    return tmp, float(n_branches * U)

branch_slide_move

branch_slide_move(tree: HifukuTree, taxon_table: TaxonTable, rng: Generator, sigma: float = 0.05) -> tuple[HifukuTree, float]

Single-branch additive slide.

Adds a Normal(0, sigma) perturbation to one randomly chosen branch length, reflecting at zero to maintain positivity. The folded-normal proposal is symmetric, so log_hastings_ratio = 0.

Parameters:

Name	Type	Description	Default
tree	HifukuTree		required
taxon_table	TaxonTable		required
rng	Generator		required
sigma	float	Standard deviation of the Normal proposal (default 0.05).	0.05

Returns:

Name	Type	Description
new_tree	HifukuTree
log_hastings_ratio	float	Always 0.0.

Notes

• Minh et al. (2020) : single-branch optimization context.

Source code in src/hifuku/moves.py

def branch_slide_move(
    tree: HifukuTree,
    taxon_table: TaxonTable,
    rng: np.random.Generator,
    sigma: float = 0.05,
) -> tuple[HifukuTree, float]:
    """
    Single-branch additive slide.

    Adds a ``Normal(0, sigma)`` perturbation to one randomly chosen branch
    length, reflecting at zero to maintain positivity. The folded-normal
    proposal is symmetric, so ``log_hastings_ratio = 0``.

    Parameters
    ----------
    tree : HifukuTree
    taxon_table : TaxonTable
    rng : np.random.Generator
    sigma : float
        Standard deviation of the Normal proposal (default 0.05).

    Returns
    -------
    new_tree : HifukuTree
    log_hastings_ratio : float
        Always 0.0.

    Notes
    -----
    - Minh et al. (2020) : single-branch optimization context.
    """
    eligible = [k for k in range(tree.n_nodes) if math.isfinite(float(tree.edge_len[k]))]
    if not eligible:
        return tree, 0.0

    k = int(rng.choice(eligible))
    delta = float(rng.normal(0.0, sigma))
    new_len = abs(float(tree.edge_len[k]) + delta)

    new_el = tree.edge_len.copy()
    new_el[k] = new_len

    tmp = HifukuTree(
        parent=tree.parent.copy(), left_child=tree.left_child.copy(),
        right_child=tree.right_child.copy(), edge_len=new_el,
        names=list(tree.names), n_nodes=tree.n_nodes,
        n_leaves=tree.n_leaves, root=tree.root,
        node_order=tree.node_order.copy() if tree.node_order is not None else None,
    )
    if tree.global_leaf_indices is not None:
        tmp.global_leaf_indices = tree.global_leaf_indices.copy()
    return tmp, 0.0