API reference

The Python API has one canonical path per name, reached through each level's package — there are no top-level re-exports. A run always starts from a simulate_* entry point and returns a *Result object that carries the true history behind the dataset.

from zombi2 import species, genomes
from zombi2.rates import scope, modifiers

result = species.simulate_species_tree(birth=1.0, death=0.3, n_extant=20, seed=1)

The reference below is generated from the source docstrings, level by level, in the same order as the user guide.

Species trees

zombi2.species.simulate_species_tree

simulate_species_tree(birth, death=0.0, *, n_extant=None, total_time=None, mass_extinctions=None, sampling=1.0, fossils=0.0, seed=None, progress=False, max_lineages=100000) -> SpeciesResult

Grow a forward birth-death tree.

birth and death are rate specs (a number, a scope wrapper, or a product with modifiers); the default scope is per lineage (each lineage speciates/dies at the base rate, so the tree grows exponentially). Yule = death=0.

Stop at exactly n_extant living lineages, or at total_time — give exactly one. n_extant is conditioned on survival: a birth-death tree can die out, so we restart (advancing the same generator) until one reaches n_extant. total_time is not conditioned. Deterministic given seed.

mass_extinctions is a list of (time, fraction_lost) pulses — e.g. [(3.0, 0.75)] culls 75% of the lineages alive at time 3.0 (time runs forward from the origin, t=0). It is a point-in-time intervention on the process (not a rate) placed on the timeline, so it needs a fixed end: give total_time (not n_extant), with each time strictly inside (0, total_time).

sampling (ρ, default 1.0) is incomplete extant sampling: each survivor is observed with probability ρ, the rest relabelled unsampled. It prunes the extant tree to the sampled survivors (the unsampled ones remain only in the complete tree). n_extant still stops at that many survivors; sampling then thins what you observe, so result.n_extant can be smaller.

fossils is a recovery rate along the branches: each branch of length L yields Poisson(fossils × L) fossils, returned as result.fossils = (lineage, time) pairs. A side output — the fossil's lineage is not removed and does not enter the extant tree.

Source code in zombi2/species/__init__.py

def simulate_species_tree(birth, death=0.0, *, n_extant=None, total_time=None,
                          mass_extinctions=None, sampling=1.0, fossils=0.0, seed=None,
                          progress=False, max_lineages=100_000) -> SpeciesResult:
    """Grow a forward birth-death tree.

    ``birth`` and ``death`` are rate specs (a number, a ``scope`` wrapper, or a product
    with modifiers); the default scope is **per lineage** (each lineage speciates/dies at
    the base rate, so the tree grows exponentially). Yule = ``death=0``.

    Stop at exactly ``n_extant`` living lineages, **or** at ``total_time`` — give exactly
    one. ``n_extant`` is **conditioned on survival**: a birth-death tree can die out, so we
    restart (advancing the same generator) until one reaches ``n_extant``. ``total_time`` is not
    conditioned. Deterministic given ``seed``.

    ``mass_extinctions`` is a list of ``(time, fraction_lost)`` pulses — e.g. ``[(3.0, 0.75)]`` culls
    75% of the lineages alive at time 3.0 (time runs forward from the origin, t=0). It is a point-in-time
    intervention on the process (not a rate) placed on the timeline, so it needs a fixed end:
    give ``total_time`` (not ``n_extant``), with each time strictly inside ``(0, total_time)``.

    ``sampling`` (ρ, default 1.0) is incomplete extant sampling: each survivor is observed with
    probability ρ, the rest relabelled ``unsampled``. It prunes the **extant tree** to the sampled
    survivors (the unsampled ones remain only in the complete tree). ``n_extant`` still stops at that
    many *survivors*; sampling then thins what you observe, so ``result.n_extant`` can be smaller.

    ``fossils`` is a recovery rate along the branches: each branch of length ``L`` yields
    ``Poisson(fossils × L)`` fossils, returned as ``result.fossils`` = ``(lineage, time)`` pairs. A
    **side output** — the fossil's lineage is not removed and does not enter the extant tree.
    """
    birth_rate = as_rate(birth, default_scope=PerLineage)
    death_rate = as_rate(death, default_scope=PerLineage)
    for label, rate in (("birth", birth_rate), ("death", death_rate)):
        # a modifier this engine does not thread would return its default factor of 1.0 — a run that
        # is quietly not the model asked for — so reject it (SPEC §5, the genome engine's discipline)
        if not isinstance(rate.scope, WIRED_SCOPES):
            raise ValueError(
                f"{label} has a {type(rate.scope).__name__} scope, but the species engine counts "
                f"lineages — use PerLineage(...) (the default, so a bare number is enough) or "
                f"Global(...) for one shared budget."
            )
        for m in rate.modifiers:
            if not isinstance(m, WIRED_MODIFIERS):
                raise ValueError(
                    f"{label} carries {type(m).__name__}, which the species engine does not "
                    f"support — OnTime (skyline), OnTotalDiversity (diversity-dependent) and "
                    f"FromParent (clade drift, ClaDS) are wired."
                )
        if _drift(rate) is not None and not isinstance(rate.scope, PerLineage):
            raise ValueError(
                f"{label} carries FromParent (per-lineage drift) but its scope is "
                f"{type(rate.scope).__name__}; a drifting rate must be per lineage — drop the "
                f"scope wrapper (per lineage is the default) or use PerLineage(...)"
            )
    if (n_extant is None) == (total_time is None):
        raise ValueError("give exactly one of n_extant or total_time")
    if n_extant is not None and (isinstance(n_extant, bool) or not isinstance(n_extant, int) or n_extant < 1):
        raise ValueError(f"n_extant must be a positive integer, got {n_extant!r}")
    if total_time is not None and (not isinstance(total_time, (int, float)) or not math.isfinite(total_time) or total_time <= 0):
        raise ValueError(f"total_time must be a positive finite number, got {total_time!r}")
    if isinstance(fossils, bool) or not isinstance(fossils, (int, float)) or not math.isfinite(fossils) or fossils < 0:
        raise ValueError(f"fossils must be a non-negative finite rate, got {fossils!r}")
    if isinstance(sampling, bool) or not isinstance(sampling, (int, float)) or not 0.0 < sampling <= 1.0:
        raise ValueError(f"sampling must be a fraction in (0, 1], got {sampling!r}")
    pulses = _mass_extinction_pulses(mass_extinctions, total_time)  # [] unless mass_extinctions given (needs total_time)

    rng = np.random.default_rng(seed)

    def _finish(tree: Tree, events: list[Event]) -> SpeciesResult:
        # observe (sampling relabels survivors) then recover fossils along the grown branches
        _apply_sampling(tree, sampling, rng)
        return SpeciesResult(tree, events, seed, _recover_fossils(tree, fossils, rng))

    if total_time is not None:
        tree, events = _grow(rng, birth_rate, death_rate, None, total_time, pulses, progress,
                             max_lineages)
        # A time-conditioned run is not conditioned on survival, so with death ≥ birth it can reach
        # total_time with nothing alive. An empty tree is not a sample anyone can use — the extant
        # tree is None and every downstream level would otherwise mistake the last-dying tip for a
        # survivor — so refuse it here rather than hand back a tree with no present.
        if not any(nd.fate == "extant" for nd in tree.nodes.values()):
            raise RuntimeError(
                f"the run went extinct before total_time={total_time:g}: no lineage is alive at the "
                f"present, so there is nothing to grow a genome, sequence or trait along. With death "
                f"close to or above birth, total extinction is likely — lower death, shorten "
                f"total_time, or use n_extant=... (which is conditioned on survival).")
        return _finish(tree, events)

    for _ in range(_MAX_ATTEMPTS):
        tree, events = _grow(rng, birth_rate, death_rate, n_extant, None, [], progress,
                             max_lineages)
        if sum(1 for nd in tree.nodes.values() if nd.fate == "extant") == n_extant:  # survivors (pre-sampling)
            return _finish(tree, events)
    raise RuntimeError(
        f"could not grow a tree to {n_extant} extant lineages in {_MAX_ATTEMPTS} attempts; "
        "birth must comfortably exceed death for large n_extant"
    )

zombi2.species.SpeciesResult dataclass

SpeciesResult(complete_tree: Tree, events: list[Event], seed: int | None, fossils: list[tuple[int, float]] = list())

What simulate_species_tree returns: the complete_tree (with the dead) and the derived extant_tree (the observed survivors), the events log (the recorded true history), the seed, and any fossils. (The record= memory dial lands with the data-heavy levels.)

n_extant property

n_extant: int

The number of observed survivors — the extant tips. Under sampling < 1 this is the sampled subset (the rest are unsampled), so it matches the extant tree's tip count.

extant_tree cached property

extant_tree: Tree | None

The survivors' tree — the complete tree pruned to extant lineages with the unifurcations suppressed (dated, bifurcating). None if nothing survived.

write

write(directory, outputs=None) -> None

Write outputs to directory, each file prefixed species_; outputs selects which (default = all applicable): "complete" → species_complete.nwk, "extant" → species_extant.nwk (if any survived), "events" → species_events.tsv (the always-recorded true history), "fossils" → species_fossils.tsv (if any recovered), "fates" → species_fates.tsv (each tip's resolved fate).

species_fates.tsv is the tip-fate table: one lineage<TAB>fate row per tip, with fate one of extant / extinct / unsampled. Fate is resolved once, at the end of the run, on the same stable n<id> that keys every other file, so it never renames anything — it is a materialised view of information the run already holds. It exists because the .nwk records only branch lengths, from which a reader cannot tell an extinct tip from a survivor that sits at the present; this table says so directly, so a downstream level can build the extant set from fate rather than guessing from tip depth.

Source code in zombi2/species/__init__.py

def write(self, directory, outputs=None) -> None:
    """Write outputs to ``directory``, each file prefixed ``species_``; ``outputs`` selects which
    (default = all applicable): ``"complete"`` → ``species_complete.nwk``, ``"extant"`` →
    ``species_extant.nwk`` (if any survived), ``"events"`` → ``species_events.tsv`` (the
    always-recorded true history), ``"fossils"`` → ``species_fossils.tsv`` (if any recovered),
    ``"fates"`` → ``species_fates.tsv`` (each tip's resolved fate).

    ``species_fates.tsv`` is the tip-fate table: one ``lineage<TAB>fate`` row per tip, with fate
    one of ``extant`` / ``extinct`` / ``unsampled``. Fate is resolved once, at the end of the run,
    on the same stable ``n<id>`` that keys every other file, so it never renames anything — it is
    a materialised view of information the run already holds. It exists because the ``.nwk`` records
    only branch lengths, from which a reader cannot tell an extinct tip from a survivor that sits at
    the present; this table says so directly, so a downstream level can build the extant set from
    fate rather than guessing from tip depth."""
    if outputs is None:
        outputs = _WRITE_OUTPUTS
    unknown = [o for o in outputs if o not in _WRITE_OUTPUTS]
    if unknown:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(_WRITE_OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    if "complete" in outputs:
        (d / "species_complete.nwk").write_text(self.complete_tree.to_newick() + "\n")
    if "extant" in outputs and self.extant_tree is not None:
        (d / "species_extant.nwk").write_text(self.extant_tree.to_newick() + "\n")
    if "events" in outputs:
        rows = ["time\tkind\tlineage\tchildren"]
        for e in self.events:
            kids = ";".join(f"n{c}" for c in e.children) if e.children else ""
            rows.append(f"{e.time:.6g}\t{e.kind}\tn{e.node}\t{kids}")
        (d / "species_events.tsv").write_text("\n".join(rows) + "\n")
    if "fossils" in outputs and self.fossils:
        rows = ["lineage\ttime"] + [f"n{i}\t{t:.6g}" for i, t in self.fossils]
        (d / "species_fossils.tsv").write_text("\n".join(rows) + "\n")
    if "fates" in outputs:
        # one row per tip (extant / extinct / unsampled); internal nodes are always speciations
        rows = ["lineage\tfate"]
        for n in sorted(self.complete_tree.leaves(), key=lambda nd: nd.id):
            rows.append(f"n{n.id}\t{n.fate}")
        (d / "species_fates.tsv").write_text("\n".join(rows) + "\n")

zombi2.species.Tree dataclass

Tree(nodes: dict[int, Node], root: int)

The complete tree: every lineage that ever lived, keyed by id, rooted at root.

leaves

leaves() -> list[Node]

Every lineage with no descendants — extant and extinct.

Source code in zombi2/tree.py

def leaves(self) -> list[Node]:
    """Every lineage with no descendants — extant **and** extinct."""
    return [n for n in self.nodes.values() if n.children is None]

extant

extant() -> list[Node]

The lineages alive at the present.

Source code in zombi2/tree.py

def extant(self) -> list[Node]:
    """The lineages alive at the present."""
    return [n for n in self.nodes.values() if n.fate == "extant"]

extinct

extinct() -> list[Node]

The lineages that died before the present.

Source code in zombi2/tree.py

def extinct(self) -> list[Node]:
    """The lineages that died before the present."""
    return [n for n in self.nodes.values() if n.fate == "extinct"]

unsampled

unsampled() -> list[Node]

Survivors not observed under incomplete sampling — kept in the complete tree (told apart by their fate) but pruned from the extant tree.

Source code in zombi2/tree.py

def unsampled(self) -> list[Node]:
    """Survivors not observed under incomplete ``sampling`` — kept in the complete tree (told
    apart by their fate) but pruned from the extant tree."""
    return [n for n in self.nodes.values() if n.fate == "unsampled"]

to_newick

to_newick() -> str

Serialise to Newick (matching tree.to_newick() elsewhere in the codebase). Each branch length is end_time - birth_time and every node — leaves and internals — is named n<id>.

The root carries a branch length like any other node: its stem, the time from the origin to the first split. A forward birth–death run starts from one lineage, so that stem is real simulated time in which events happen, and writing )n0; would silently discard it — for a tree whose crown comes late, a large fraction of its history. It is emitted as )n0:<stem>; and :func:read_newick reads it back.

Source code in zombi2/tree.py

def to_newick(self) -> str:
    """Serialise to Newick (matching ``tree.to_newick()`` elsewhere in the codebase). Each
    branch length is ``end_time - birth_time`` and every node — leaves and internals — is named
    ``n<id>``.

    The root carries a branch length like any other node: its **stem**, the time from the origin
    to the first split. A forward birth–death run starts from one lineage, so that stem is real
    simulated time in which events happen, and writing ``)n0;`` would silently discard it — for a
    tree whose crown comes late, a large fraction of its history. It is emitted as ``)n0:<stem>;``
    and :func:`read_newick` reads it back."""

    def emit(i: int) -> str:
        node = self.nodes[i]
        bl = node.end_time - node.birth_time
        if node.children is None:
            return f"n{i}:{bl:.6g}"
        inner = ",".join(emit(c) for c in node.children)
        return f"({inner})n{i}:{bl:.6g}"

    root = self.nodes[self.root]
    stem = root.end_time - root.birth_time
    if root.children is None:
        return f"n{self.root}:{stem:.6g};"
    return f"({','.join(emit(c) for c in root.children)})n{self.root}:{stem:.6g};"

zombi2.species.Node dataclass

Node(id: int, parent: int | None, birth_time: float, end_time: float | None = None, children: tuple[int, int] | None = None, fate: str = 'alive')

One lineage segment: born at birth_time, ended at end_time by a split, a death, or reaching the present. A split has two children; a leaf has none.

zombi2.species.Event dataclass

Event(time: float, kind: str, node: int, children: tuple[int, int] | None = None)

A recorded event in the true history: a speciation (with its two children) or an extinction.

zombi2.species.prune

prune(tree: Tree, keep: str = 'extant') -> Tree | None

Prune the complete tree to a kept set (matching prune(tree, keep=...) in the codebase): drop the pruned subtrees and suppress the unifurcations they leave behind, giving a dated, bifurcating tree. Branch lengths merge across suppressed nodes; None if nothing is kept.

keep="extant" (default) keeps the survivors — the extant tree. ("sampled", the fossil/serially-sampled tree, arrives with the sampling and fossils slices; it raises for now.)

Source code in zombi2/tree.py

def prune(tree: Tree, keep: str = "extant") -> Tree | None:
    """Prune the complete tree to a kept set (matching ``prune(tree, keep=...)`` in the codebase):
    drop the pruned subtrees and suppress the unifurcations they leave behind, giving a dated,
    bifurcating tree. Branch lengths merge across suppressed nodes; ``None`` if nothing is kept.

    ``keep="extant"`` (default) keeps the survivors — the extant tree. (``"sampled"``, the
    fossil/serially-sampled tree, arrives with the sampling and fossils slices; it raises for now.)"""
    if keep != "extant":
        raise ValueError(f"keep must be 'extant' until the sampling/fossils slices land, got {keep!r}")
    nodes = tree.nodes
    surviving: dict[int, bool] = {}
    for i in sorted(nodes, reverse=True):  # children have higher ids → processed before parents
        nd = nodes[i]
        surviving[i] = nd.fate == "extant" if nd.children is None else any(surviving[c] for c in nd.children)
    if not any(surviving.values()):
        return None

    def surv_children(i: int) -> list[int]:
        nd = nodes[i]
        return [] if nd.children is None else [c for c in nd.children if surviving[c]]

    # keep the extant leaves and the genuine bifurcations (≥2 surviving children)
    kept = {i for i in nodes
            if (nodes[i].children is None and nodes[i].fate == "extant") or len(surv_children(i)) >= 2}

    new: dict[int, Node] = {}
    ext_root: int | None = None
    for i in kept:
        p = nodes[i].parent  # walk up to the nearest kept ancestor
        while p is not None and p not in kept:
            p = nodes[p].parent
        branch_start = nodes[p].end_time if p is not None else 0.0  # merge the suppressed edges
        new[i] = Node(i, p, branch_start, nodes[i].end_time, None, nodes[i].fate)
        if p is None:
            ext_root = i
    for i in sorted(kept):  # rebuild children from parents, in id order for a stable Newick
        p = new[i].parent
        if p is not None:
            existing = new[p].children
            new[p].children = (i,) if existing is None else existing + (i,)

    return Tree(new, ext_root)

Genomes

The genome level has three resolutions — family ⊂ ordered ⊂ nucleotide — one entry point each.

zombi2.genomes.simulate_genomes_family

simulate_genomes_family(tree, *, duplication=0.0, transfer=0.0, loss=0.0, origination=0.0, transfer_to='uniform', replacement=False, self_transfer=False, initial_families=100, family_names=None, family_speed=None, max_family_size=10.0, seed=None, parallel=False, stream_to=None, outputs=None, progress=False) -> 'FamilyGenomesResult | StreamedRun'

Evolve a multiset of gene families along a species tree by duplication, transfer, loss, and origination.

tree is the complete species tree (a :class:~zombi2.species_tree.Tree, or a :class:~zombi2.species_tree.SpeciesResult whose complete_tree is used). Genomes evolve on every lineage, extant and extinct alike, so the true gene-tree history is complete and a transfer can arrive "from the dead"; the observed genomes are the extant tips.

Rates (each a scope(base) × modifiers spec): duplication/transfer/loss default per copy, origination per lineage. When a transfer fires it moves a copy from a uniformly-chosen donor copy to a recipient lineage alive at that instant, chosen by transfer_to — "uniform" (any other contemporaneous lineage), "distance" / Distance(decay=) (closer relatives likelier), Clades({...}, Between({...})) (weighted by the donor's and recipient's named clade, so transfer can run between two clades — see below), or mod.DrivenBy(source, mapping) (weighted by another level; see below). replacement=True overwrites a homologous copy in the recipient (additive fallback if it has none); self_transfer=True lets a lineage donate to itself. The root starts with initial_families families of one copy each, recorded as originations at the crown. family_names=["toxin", …] additionally declares named families — each gets a normal (integer) family id, but its name is remembered in result.family_names so you can track a specific family (result.has_family(node, "toxin")); this is the handle a joint DrivenBy("genomes:toxin", …) reads. Deterministic given seed.

Conditioning (a trait drives a rate). Any of the four rates may be driven by another level — loss = 0.25 * mod.DrivenBy("habitat.tsv", {"aquatic": 3.0, "terrestrial": 1.0}) scales each lineage's loss by the habitat on that branch, read from a driver file grown first (traits.simulate_discrete(...).write(dir, outputs=("driver",))). A driven rate is then per-lineage: it is summed over the living lineages (each with its own copy count and driver value), the affected lineage is drawn weighted by its rate, and the Gillespie steps at every mid-branch switch of the driver (SPEC §2). For transfer the affected lineage is the donor, so a driven transfer says how often a lineage donates.

Conditioning (a trait drives who receives). transfer_to = mod.DrivenBy(source, mapping) is the other half, and a different model: the mapping's numbers are per-candidate weights, not rate multipliers, so they leave the total amount of transfer alone and only redistribute it (SPEC §5, the choice slot). Candidate lineage k gets weight mapping(driver value on k now) and receives with probability w_k / Σw — five candidates at weight 1 and five at weight 2 send two thirds of transfers to the weight-2 group. Weight 0 means "cannot receive"; when every candidate weighs 0 the transfer does not happen (see :func:_do_transfer). The two couplings are independent and may be used together or apart.

parallel opts into a separate engine that evolves the families concurrently, one per worker process — the families are independent (a transfer roams a copy across lineages but never mixes two families), so it enumerates every family's origination first and then evolves each on its own. It is worker-count invariant (each family draws from its own stream spawned off seed) but gives a different-though-equally-valid draw than the serial default for a given seed. False (default) runs the serial loop; True uses every core; an int sets the worker count. A driven rate or transfer_to (DrivenBy) has no per-family engine yet — the run announces this and falls back to the serial loop. The gain is real but modest (a merge over the whole event log stays serial), so a handful of workers is the sweet spot; unlike the sequences level it does not scale far. Because it spawns processes, a calling script must guard its entry with if __name__ == "__main__": (the zombi2 CLI already does).

stream_to=DIR takes the same engine to the many-families regime: each family is written straight to disk as it finishes — no whole-run merge, no run held in memory (a run that fills gigabytes in memory streams in tens of megabytes) — and a light :class:~zombi2.genomes.StreamedRun handle comes back instead of a FamilyGenomesResult. outputs= picks which files, exactly as :meth:FamilyGenomesResult.write takes them (default: all of them). It is the per-family engine, so a driven rate raises here rather than falling back (that would defeat the point), and outputs without stream_to is an error.

Source code in zombi2/genomes/family.py

@without_cyclic_gc
def simulate_genomes_family(tree, *, duplication=0.0, transfer=0.0, loss=0.0, origination=0.0,
                            transfer_to="uniform", replacement=False, self_transfer=False,
                            initial_families=100, family_names=None, family_speed=None,
                            max_family_size=10.0, seed=None, parallel=False, stream_to=None,
                            outputs=None, progress=False) -> "FamilyGenomesResult | StreamedRun":
    """Evolve a multiset of gene families along a species tree by duplication, transfer, loss, and
    origination.

    ``tree`` is the **complete** species tree (a :class:`~zombi2.species_tree.Tree`, or a
    :class:`~zombi2.species_tree.SpeciesResult` whose ``complete_tree`` is used). Genomes evolve on
    **every** lineage, extant and extinct alike, so the true gene-tree history is complete and a
    transfer can arrive "from the dead"; the observed genomes are the extant tips.

    Rates (each a ``scope(base) × modifiers`` spec): ``duplication``/``transfer``/``loss`` default
    **per copy**, ``origination`` **per lineage**. When a transfer fires it moves a copy from a
    uniformly-chosen donor copy to a recipient lineage alive at that instant, chosen by
    ``transfer_to`` — ``"uniform"`` (any other contemporaneous lineage), ``"distance"`` /
    ``Distance(decay=)`` (closer relatives likelier), ``Clades({...}, Between({...}))`` (weighted by
    the donor's and recipient's **named clade**, so transfer can run *between* two clades — see below),
    or ``mod.DrivenBy(source, mapping)`` (weighted by another level; see below). ``replacement=True``
    overwrites a homologous
    copy in the recipient (additive fallback if it has none); ``self_transfer=True`` lets a lineage
    donate to itself. The root starts with ``initial_families`` families of one copy each, recorded
    as originations at the crown. ``family_names=["toxin", …]`` additionally declares **named** families —
    each gets a normal (integer) family id, but its name is remembered in ``result.family_names`` so
    you can track a specific family (``result.has_family(node, "toxin")``); this is the handle a joint
    ``DrivenBy("genomes:toxin", …)`` reads. Deterministic given ``seed``.

    **Conditioning (a trait drives a rate).** Any of the four rates may be *driven by another level* —
    ``loss = 0.25 * mod.DrivenBy("habitat.tsv", {"aquatic": 3.0, "terrestrial": 1.0})`` scales each
    lineage's loss by the habitat on that branch, read from a driver file grown first
    (``traits.simulate_discrete(...).write(dir, outputs=("driver",))``). A driven rate is then
    *per-lineage*: it is summed over the living lineages (each with its own copy count and driver
    value), the affected lineage is drawn weighted by its rate, and the Gillespie steps at every
    mid-branch switch of the driver (SPEC §2). For ``transfer`` the affected
    lineage is the **donor**, so a driven ``transfer`` says how often a lineage *donates*.

    **Conditioning (a trait drives who receives).** ``transfer_to = mod.DrivenBy(source, mapping)`` is
    the other half, and a different model: the mapping's numbers are per-candidate **weights**, not
    rate multipliers, so they leave the total amount of transfer alone and only redistribute it
    (SPEC §5, the choice slot). Candidate lineage ``k`` gets weight ``mapping(driver value on k now)``
    and receives with probability ``w_k / Σw`` — five candidates at weight 1 and five at weight 2 send
    two thirds of transfers to the weight-2 group. Weight 0 means "cannot receive"; when every
    candidate weighs 0 the transfer does not happen (see :func:`_do_transfer`). The two couplings are
    independent and may be used together or apart.

    ``parallel`` opts into a **separate** engine that evolves the families concurrently, one per worker
    process — the families are independent (a transfer roams a copy across lineages but never mixes two
    families), so it enumerates every family's origination first and then evolves each on its own. It is
    worker-count invariant (each family draws from its own stream spawned off ``seed``) but gives a
    different-though-equally-valid draw than the serial default for a given seed. ``False`` (default)
    runs the serial loop; ``True`` uses every core; an ``int`` sets the worker count. A **driven** rate
    or ``transfer_to`` (``DrivenBy``) has no per-family engine yet — the run announces this and falls
    back to the serial loop. The gain is real but modest (a merge over the whole event log stays
    serial), so a handful of workers is the sweet spot; unlike the sequences level it does not scale far.
    Because it spawns processes, a calling script must guard its entry with ``if __name__ ==
    "__main__":`` (the ``zombi2`` CLI already does).

    ``stream_to=DIR`` takes the same engine to the many-families regime: each family is written straight
    to disk as it finishes — no whole-run merge, no run held in memory (a run that fills gigabytes in
    memory streams in tens of megabytes) — and a light :class:`~zombi2.genomes.StreamedRun` handle comes
    back instead of a ``FamilyGenomesResult``. ``outputs=`` picks which files, exactly as
    :meth:`FamilyGenomesResult.write` takes them (default: all of them). It is the per-family engine, so a
    driven rate **raises** here rather than falling back (that would defeat the point), and ``outputs``
    without ``stream_to`` is an error.
    """
    tree = tree.complete_tree if isinstance(tree, SpeciesResult) else tree
    dup = as_rate(duplication, default_scope=PerCopy)
    tra = as_rate(transfer, default_scope=PerCopy)
    los = as_rate(loss, default_scope=PerCopy)
    org = as_rate(origination, default_scope=PerLineage)
    # this slice wires only the default scope (D/T/L per copy, origination per lineage). The wired
    # modifiers are OnTime (skyline) and DrivenBy (a conditioned/joint driver). A non-default scope
    # would set the *total* rate one way while the engine still picks the affected copy/lineage the
    # default way — a silent mismatch (e.g. a PerCopy origination is base×0 copies, a no-op) — so
    # reject it. DrivenBy is a per-lineage driver, wired for all four events; on transfer the driven
    # lineage is the DONOR (who receives is the separate transfer_to slot, below).
    for label, rate, want in (("duplication", dup, PerCopy), ("transfer", tra, PerCopy),
                              ("loss", los, PerCopy), ("origination", org, PerLineage)):
        if not isinstance(rate.scope, want):
            raise ValueError(
                f"{label} has a {type(rate.scope).__name__} scope, but the family genome engine "
                f"wires only {want.__name__} for {label} this slice — scope overrides are a later slice."
            )
        for m in rate.modifiers:
            if isinstance(m, ByFamily) and label == "origination":
                raise ValueError(
                    "origination carries ByFamily, but origination is the rate at which families are "
                    "CREATED — when it is read there is no family yet to have drawn a factor for. "
                    "Put ByFamily on duplication, transfer or loss, or use family_speed= for a "
                    "family-wide tempo.")
            if isinstance(m, DrivenBy) and isinstance(m.mapping, Between):
                raise ValueError(
                    f"{label} carries DrivenBy(…, Between(…)); a Between kernel is donor-conditioned — "
                    f"it weights a recipient by the (donor, recipient) group pair — so it belongs only "
                    f"in transfer_to (who RECEIVES), never in a rate (which has no donor to condition "
                    f"on). Drive a rate with a Table (a plain dict).")
            if isinstance(m, (OnTime, DrivenBy, ByFamily)):
                continue
            raise ValueError(
                f"{label} carries {type(m).__name__}, which the family genome engine does not "
                f"support yet — only OnTime (skyline), DrivenBy (a conditioned/joint driver) and "
                f"ByFamily (per-family heterogeneity) are wired. Clade drift is a later slice."
            )
    if any(isinstance(m, ByFamily) for rate in (dup, tra, los) for m in rate.modifiers) and \
            any(isinstance(m, DrivenBy) for rate in (dup, tra, los, org) for m in rate.modifiers):
        raise ValueError(
            "ByFamily and DrivenBy on the same run is a later slice: one weights lineages by a "
            "driver and the other weights copies by their family, and combining them means "
            "weighting by the product. Use one or the other for now.")
    if family_speed is not None and not isinstance(family_speed, ByFamily):
        raise ValueError(
            f"family_speed takes a ByFamily modifier — family_speed=mod.ByFamily(spread=0.5) — "
            f"got {family_speed!r}. It is the family-wide slot: one draw per family scaling every "
            f"rate that family has.")
    if transfer_to == "distance":
        transfer_to = Distance()
    if isinstance(transfer_to, Rate):
        # `1.0 * mod.DrivenBy(...)` — the rate spelling, in a slot that is not a rate. Say so rather
        # than let a Rate fall through to the generic "must be …" message.
        raise ValueError(
            "transfer_to takes the DrivenBy modifier on its own, not a rate — write "
            "transfer_to=mod.DrivenBy(source, {...}) with no base number. In this slot the mapping's "
            "numbers are relative WEIGHTS over the candidate recipients (normalised), not a rate "
            "multiplier: they change who receives, never how much transfer happens."
        )
    if isinstance(transfer_to, (list, tuple)):
        raise ValueError(
            "transfer_to takes one recipient rule, not several — combining Distance (relatedness) "
            "with a DrivenBy weighting is a later slice. Give 'uniform', 'distance' / "
            "Distance(decay=), or mod.DrivenBy(source, {...})."
        )
    if transfer_to != "uniform" and not isinstance(transfer_to, (Distance, DrivenBy, Clades)):
        raise ValueError(
            f"transfer_to must be 'uniform', 'distance' / Distance(decay=), "
            f"Clades({{...}}, Between({{...}})) (weight by named clade), or "
            f"mod.DrivenBy(source, {{...}}) (a recipient weight driven by another level), "
            f"got {transfer_to!r}")
    if isinstance(initial_families, bool) or not isinstance(initial_families, int) or initial_families < 0:
        raise ValueError(f"initial_families must be a non-negative integer, got {initial_families!r}")
    family_names = list(family_names) if family_names is not None else []
    for name in family_names:
        if not isinstance(name, str) or not name.strip():
            raise ValueError(f"family_names must be a list of non-empty family names (strings), got {name!r}")
    if len(set(family_names)) != len(family_names):
        raise ValueError(f"family names must be unique, got {family_names}")

    # A family's copies in one genome are capped. Growth compounds — a duplication rate above the
    # loss rate multiplies without bound — so a run needs a ceiling somewhere. An int is that number
    # of copies; a float (the default, 10.0) is that multiple of the lineages in the complete tree,
    # so the bound travels with the size of the run. Refusing an event on a condition that depends
    # only on the current state is Poisson thinning, so what is kept is exactly the process whose
    # duplication rate is zero for a family already at its quota — a declared ceiling, not a
    # truncated run. ``None`` removes it.
    cap = resolve_max_family_size(max_family_size, len(tree.nodes))

    # Parallel is a *separate* engine (opt-in): families are independent, so it evolves them one per
    # process (SPEC-style — serial by default). `stream_to` takes the same engine one step further —
    # each family is written straight to disk and a light StreamedRun handle comes back, so a run of a
    # million families never has to fit in memory (`outputs` picks which files, as `.write` does). Both
    # handle everything but a driven rate / recipient; there the parallel path prints why and returns
    # None so this serial reference loop runs unchanged (decision A), while a streamed run raises (it
    # cannot fall back without pulling the whole thing into memory). Everything above is shared
    # validation, so bad input still raises the same way.
    if outputs is not None and stream_to is None:
        raise ValueError(
            "outputs applies to a streamed run (stream_to=DIR), which writes the files itself; for an "
            "in-memory run choose them when you call result.write(outputs=...).")
    if parallel or stream_to is not None:
        from ._perfamily import run_parallel_family
        result = run_parallel_family(
            tree, dup=dup, tra=tra, los=los, org=org, transfer_to=transfer_to,
            replacement=replacement, self_transfer=self_transfer, initial_families=initial_families,
            family_names=family_names, family_speed=family_speed, cap=cap, seed=seed, parallel=parallel,
            progress=progress, stream_to=stream_to, outputs=outputs)
        if result is not None:
            return result

    rng = np.random.default_rng(seed)
    copy_counter = 0
    family_counter = 0

    def new_copy(family: int) -> GeneCopy:
        nonlocal copy_counter
        c = GeneCopy(copy_counter, family)
        copy_counter += 1
        return c

    # Per-family multipliers, drawn once when a family is minted and then fixed for its whole life:
    # family_speed scales every rate that family has (one draw), and a ByFamily on a single rate
    # varies that rate on its own (a separate draw). Placement is what decides whether a family's
    # rates move together. Empty unless one of them is used, and then the engine takes its weighted
    # path; a run without either draws nothing here and is byte-identical to before.
    fam_by = {"duplication": next((m for m in dup.modifiers if isinstance(m, ByFamily)), None),
              "transfer": next((m for m in tra.modifiers if isinstance(m, ByFamily)), None),
              "loss": next((m for m in los.modifiers if isinstance(m, ByFamily)), None)}
    any_family = family_speed is not None or any(fam_by.values())
    fam_mult: dict[str, dict[int, float]] = {key: {} for key in fam_by}

    def new_family() -> int:
        nonlocal family_counter
        f = family_counter
        family_counter += 1
        if any_family:
            speed = family_speed.draw(rng) if family_speed is not None else 1.0
            for key, m in fam_by.items():
                fam_mult[key][f] = speed * (m.draw(rng) if m is not None else 1.0)
        return f

    depth = mean_root_to_tip(tree)  # timescale for Distance weighting (unused by "uniform")
    # a Clades rule paints every lineage with its clade once (membership is a fact about the tree, not
    # a driver, so it is constant along a branch and adds no Gillespie breakpoints). A kernel naming
    # only absent groups weights every candidate at its default — secretly uniform — so refuse it here.
    group_of = resolve_groups(tree, transfer_to.groups) if isinstance(transfer_to, Clades) else None
    if group_of is not None:
        check_kernel_fires(transfer_to.between, set(group_of.values()), source_label="clades")
    schedule = sorted((tree.nodes[i].end_time, i) for i in tree.nodes)  # (end_time, node_id)

    root = tree.nodes[tree.root]
    t = root.birth_time
    alive: list[int] = []
    gen: list[list[GeneCopy]] = []
    pos: dict[int, int] = {}
    genomes: dict[int, tuple[GeneCopy, ...]] = {}
    events: list[Event] = []
    enter(alive, gen, pos, root.id, [])
    for _ in range(initial_families):  # lay down the crown as originations at t = root.birth_time
        _originate(gen[0], root, t, events, new_copy, new_family)
    named: dict[str, int] = {}  # a minted id per declared name (so GeneCopy.family stays an int)
    for name in family_names:
        fid = new_family()
        named[name] = fid
        c = new_copy(fid)
        gen[0].append(c)
        events.append(Event(t, "origination", root.id, fid, c.id))
    total_copies = len(gen[0])
    initial_genome = tuple(gen[0])   # the run's starting genome: a snapshot before the stem runs

    # conditioning: a rate carrying DrivenBy reads a driver per lineage. Resolve each driver once into
    # a DriverTrajectory (value + next-switch lookups, keyed by the shared species node id) — from a
    # file (str source) or an in-memory trait result (object source). No driven rate ⇒ this is empty
    # and the loop stays byte-identical to an uncoupled run.
    dup_mods, los_mods = _driven_mods(dup), _driven_mods(los)
    org_mods, tra_mods = _driven_mods(org), _driven_mods(tra)
    by_key = {}  # driver key → its source (deduped, so a driver shared across rates resolves once)
    for m in (*dup_mods, *los_mods, *org_mods, *tra_mods):
        by_key.setdefault(m.key, m.source)
    rate_keys = list(by_key)     # the drivers that move a RATE: they alone set the Gillespie horizon
    to_mod = transfer_to if isinstance(transfer_to, DrivenBy) else None
    if to_mod is not None:       # the transfer_to driver is read only at the instant a transfer fires
        by_key.setdefault(to_mod.key, to_mod.source)
    resolved = {}
    if by_key:
        from ..rates.driver import check_mapping_fires, resolve_driver
        resolved = {key: resolve_driver(src, tree) for key, src in by_key.items()}
        # a mapping whose states never occur in the driver leaves every lineage at the default factor,
        # so the rate is never driven and the run is secretly the uncoupled model — refuse it here,
        # naming the driver, rather than let it pass as a coupled run
        for m in (*dup_mods, *los_mods, *org_mods, *tra_mods, *( (to_mod,) if to_mod else () )):
            label = m.source if isinstance(m.source, str) else f"<{type(m.source).__name__}>"
            check_mapping_fires(m.mapping, resolved[m.key].states(), source_label=label)
    # a driven transfer_to changes no rate — the weights are evaluated when a transfer fires, so the
    # recipient driver deliberately stays OUT of trajs (no per-lineage rate weights, no extra horizon
    # breakpoints for it). Only the rate drivers make the loop per-lineage.
    trajs = {key: resolved[key] for key in rate_keys}
    to_traj = resolved[to_mod.key] if to_mod is not None else None
    if to_mod is not None and isinstance(to_mod.mapping, Between):  # a donor-conditioned trait kernel:
        label = to_mod.source if isinstance(to_mod.source, str) else f"<{type(to_mod.source).__name__}>"
        check_kernel_fires(to_mod.mapping, to_traj.states(), source_label=label)
    any_driven = bool(rate_keys)

    # the species tree's schedule is the run's spine: one entry per speciation/extinction, so how
    # far through it we are is how far through the tree the genomes have got
    bar = progress_bar(len(schedule), "genomes", unit="branch", enabled=progress)
    si = 0
    while si < len(schedule):
        bar.to(si)
        n = total_copies
        k_alive = len(alive)
        ctx = {"copies": n, "lineages": k_alive, "time": t}
        can_xfer = n > 0 and (k_alive >= 2 or self_transfer)  # a recipient must be able to exist
        # a driven rate is per-lineage: sum its effective rate over the living lineages (each read with
        # its own copy count and its branch's driver value), keeping the weights for the affected-lineage
        # pick — the species_tree._grow shape. An undriven rate stays pooled (one .effective, uniform
        # pick), so a run with no coupling is byte-identical to before. For transfer the affected
        # lineage is the donor, so a driven transfer weights who donates.
        w_dup = w_los = w_org = w_tra = None
        if any_family:
            # A per-copy rate pools over copies, so with per-family multipliers the total is the
            # unit rate times the sum of those multipliers over the live copies — and the copy has
            # to be drawn with the same weights, or the rates would say one thing and the picking
            # another. Summed per lineage, so the existing weighted-lineage pick can be reused.
            fw = {key: [sum(fam_mult[key][c.family] for c in gen[k]) for k in range(k_alive)]
                  for key in fam_mult}
            unit = {"duplication": dup.effective(copies=1, lineages=1, time=t),
                    "loss": los.effective(copies=1, lineages=1, time=t),
                    "transfer": tra.effective(copies=1, lineages=1, time=t) if can_xfer else 0.0}
            w_dup = [unit["duplication"] * s for s in fw["duplication"]]
            w_los = [unit["loss"] * s for s in fw["loss"]]
            if can_xfer:
                w_tra = [unit["transfer"] * s for s in fw["transfer"]]
        if any_driven:
            drivers = [{key: trajs[key].value(alive[k], t) for key in trajs} for k in range(k_alive)]
            if dup_mods:
                w_dup = [dup.effective(copies=len(gen[k]), lineages=1, time=t, drivers=drivers[k])
                         for k in range(k_alive)]
            if los_mods:
                w_los = [los.effective(copies=len(gen[k]), lineages=1, time=t, drivers=drivers[k])
                         for k in range(k_alive)]
            if org_mods:
                w_org = [org.effective(copies=len(gen[k]), lineages=1, time=t, drivers=drivers[k])
                         for k in range(k_alive)]
            if tra_mods and can_xfer:
                w_tra = [tra.effective(copies=len(gen[k]), lineages=1, time=t, drivers=drivers[k])
                         for k in range(k_alive)]
        r_dup = sum(w_dup) if w_dup is not None else (dup.effective(**ctx) if n else 0.0)
        r_los = sum(w_los) if w_los is not None else (los.effective(**ctx) if n else 0.0)
        r_org = sum(w_org) if w_org is not None else org.effective(**ctx)
        r_tra = sum(w_tra) if w_tra is not None else (tra.effective(**ctx) if can_xfer else 0.0)
        total = r_dup + r_los + r_org + r_tra

        next_species = schedule[si][0]  # the tree's own next event: who is alive changes only here
        horizon = min(next_species, dup.next_change(t), los.next_change(t),
                      org.next_change(t), tra.next_change(t))
        if any_driven:  # a driven rate also changes when the driver switches mid-branch — step there
            driver_next = min((trajs[key].next_change(alive[k], t) for key in trajs
                               for k in range(k_alive)), default=math.inf)
            horizon = min(horizon, driver_next)

        if total > 0.0:
            t_ev = t + float(rng.exponential(1.0 / total))
            if t_ev < horizon:  # a genome event fires before the alive set or the rate changes
                t = t_ev
                r = float(rng.random()) * total
                if r < r_dup:
                    if w_dup is not None:  # weighted lineage, then a copy within it
                        k = weighted_index(rng, w_dup, r_dup)
                        j = (_pick_copy_by_family(rng, gen[k], fam_mult["duplication"])
                             if any_family else int(rng.integers(len(gen[k]))))
                    else:
                        k, j = _pick_copy(rng, gen, n)
                    if not _at_cap(gen[k], gen[k][j].family, cap):
                        _duplicate(gen[k], j, tree.nodes[alive[k]], t, events, new_copy)
                        total_copies += 1
                elif r < r_dup + r_los:
                    if w_los is not None:
                        k = weighted_index(rng, w_los, r_los)
                        j = (_pick_copy_by_family(rng, gen[k], fam_mult["loss"])
                             if any_family else int(rng.integers(len(gen[k]))))
                    else:
                        k, j = _pick_copy(rng, gen, n)
                    _lose_at(gen[k], j, tree.nodes[alive[k]], t, events)
                    total_copies -= 1
                elif r < r_dup + r_los + r_org:
                    k = (weighted_index(rng, w_org, r_org) if w_org is not None
                         else int(rng.integers(k_alive)))  # origination is per lineage
                    _originate(gen[k], tree.nodes[alive[k]], t, events, new_copy, new_family)
                    total_copies += 1
                else:
                    if w_tra is not None:  # driven: weighted DONOR lineage, then a uniform copy in it
                        kd = weighted_index(rng, w_tra, r_tra)
                        if not gen[kd]:    # only via _weighted_index's r == total float guard: a
                            # zero-weight lineage has no copies to donate, so take the heaviest instead
                            kd = max(range(k_alive), key=lambda k: w_tra[k])
                        jd = (_pick_copy_by_family(rng, gen[kd], fam_mult["transfer"])
                              if any_family else int(rng.integers(len(gen[kd]))))
                    else:
                        kd, jd = _pick_copy(rng, gen, n)
                    total_copies += _do_transfer(rng, tree, alive, gen, kd, jd, t, events, new_copy,
                                                 transfer_to, replacement, self_transfer, depth,
                                                 to_traj, cap, group_of)
                continue

        if horizon == next_species:  # advance to the tree's next event(s); process the whole tie-batch
            t = next_species
            while si < len(schedule) and schedule[si][0] == t:
                i = schedule[si][1]
                g = gen[pos[i]]
                genomes[i] = tuple(g)  # finalise this lineage (extant, extinct, or unsampled)
                total_copies -= len(g)
                retire(alive, gen, pos, pos[i])
                node = tree.nodes[i]
                if node.children is not None:  # a speciation: each gene re-ids into each daughter
                    for c in node.children:
                        child_genome = []
                        for old in g:  # ZOMBI1: the gene ends here and continues under a fresh id
                            nc = new_copy(old.family)
                            child_genome.append(nc)
                            events.append(Event(t, "speciation", c, old.family, nc.id, parent=old.id))
                        enter(alive, gen, pos, c, child_genome)
                        total_copies += len(child_genome)
                si += 1
        else:
            t = horizon  # a skyline breakpoint: advance and re-evaluate the (now changed) rate

    bar.close()
    return FamilyGenomesResult(tree, genomes, events, seed, named, initial_genome)

zombi2.genomes.simulate_genomes_ordered

simulate_genomes_ordered(tree, *, duplication=0.0, transfer=0.0, loss=0.0, origination=0.0, inversion=0.0, transposition=0.0, translocation=0.0, chromosomes=1, topology='circular', fission=0.0, fusion=0.0, chromosome_origination=0.0, chromosome_loss=0.0, duplication_extent=None, loss_extent=None, transfer_extent=None, inversion_extent=None, transposition_extent=None, translocation_extent=None, inversion_probability=0.0, transfer_to='uniform', replacement=False, self_transfer=False, initial_families=100, family_names=None, family_speed=None, max_family_size=10.0, seed=None, progress=False) -> OrderedGenomesResult

Evolve ordered genomes — genes with a position and an orientation, on chromosomes — along a species tree, by the D/T/L/O core plus segmental rearrangements and the chromosome tier.

Every gene-level event acts on an extent — a run of consecutive genes (the ZOMBI1 model): duplication copies the run in tandem, loss removes it, transfer sends it to a contemporaneous recipient as a block, inversion reverses it (flipping strands), transposition relocates it elsewhere on the same chromosome, and translocation moves it to a different chromosome. The run's length is drawn per event from <event>_extent (a distribution, default Geometric(mean=1) — usually a single gene; dial the mean up for larger blocks). origination is the exception: a family is born once, a single gene, no extent. transposition and translocation land the moved block inverted with probability inversion_probability.

Where a run stops is set by the chromosome's topology. A run goes rightwards from the gene it starts at. On a "circular" chromosome there are no ends, so a run that reaches the last gene continues from the first, and only the whole chromosome bounds it; on a "linear" one the run stops at the last gene. So on a circular chromosome every gene is covered by segmental events at the same rate, and the nominal mean extent is the realised one.

Scopes follow the cross-level grammar, which counts an event per the thing it acts on: the gene-level events — duplication/transfer/loss and the rearrangements inversion/transposition/translocation — are per copy, since each acts on a run of genes that starts at one of them; the chromosome tier fission/fusion/chromosome_loss is per chromosome; and the two events that make something from nothing, origination/chromosome_origination, are per lineage. The run starts with chromosomes chromosomes of the given topology, across which the initial_families founding genes are dealt round-robin; family_names=["toxin", …] additionally declares named families (remembered in result.family_names for result.has_family(node, "toxin")), as in the family core; transfer_to / replacement / self_transfer behave as in the family core.

The chromosome tier changes chromosome number: fission (split), fusion (merge — the reticulation), chromosome_origination (a de-novo replicon), chromosome_loss (a whole chromosome and its genes die; never the genome's last). Chromosomes carry identity — re-minted at every event that reshapes them — so chromosome_events is the true reticulating chromosome genealogy, rooted at the initial and de-novo originations. Deterministic given seed.

Source code in zombi2/genomes/ordered.py

@without_cyclic_gc
def simulate_genomes_ordered(tree, *, duplication=0.0, transfer=0.0, loss=0.0, origination=0.0,
                             inversion=0.0, transposition=0.0, translocation=0.0,
                             chromosomes=1, topology="circular",
                             fission=0.0, fusion=0.0, chromosome_origination=0.0, chromosome_loss=0.0,
                             duplication_extent=None, loss_extent=None, transfer_extent=None,
                             inversion_extent=None, transposition_extent=None,
                             translocation_extent=None, inversion_probability=0.0,
                             transfer_to="uniform", replacement=False, self_transfer=False,
                             initial_families=100, family_names=None, family_speed=None,
                             max_family_size=10.0, seed=None,
                             progress=False) -> OrderedGenomesResult:
    """Evolve ordered genomes — genes with a position and an orientation, on chromosomes — along a
    species tree, by the D/T/L/O core plus segmental rearrangements and the chromosome tier.

    **Every gene-level event acts on an *extent*** — a run of consecutive genes (the ZOMBI1 model):
    ``duplication`` copies the run in tandem, ``loss`` removes it, ``transfer`` sends it to a
    contemporaneous recipient as a block, ``inversion`` reverses it (flipping strands), ``transposition``
    relocates it elsewhere on the same chromosome, and ``translocation`` moves it to a different
    chromosome. The run's length is drawn per event from ``<event>_extent`` (a distribution,
    default ``Geometric(mean=1)`` — usually a single gene; dial the mean up for larger blocks).
    ``origination`` is the exception: a family is born once, a single gene, no extent.
    ``transposition`` and ``translocation`` land the moved block inverted with probability
    ``inversion_probability``.

    **Where a run stops is set by the chromosome's ``topology``.** A run goes rightwards from the gene
    it starts at. On a ``"circular"`` chromosome there are no ends, so a run that reaches the last
    gene continues from the first, and only the whole chromosome bounds it; on a ``"linear"`` one the
    run stops at the last gene. So on a circular chromosome every gene is covered by segmental events
    at the same rate, and the nominal mean extent is the realised one.

    Scopes follow the cross-level grammar, which counts an event per the thing it acts on: the
    gene-level events — ``duplication``/``transfer``/``loss`` and the rearrangements
    ``inversion``/``transposition``/``translocation`` — are **per copy**, since each acts on a run of
    genes that starts at one of them; the chromosome tier ``fission``/``fusion``/``chromosome_loss``
    is **per chromosome**; and the two events that make something from nothing,
    ``origination``/``chromosome_origination``, are **per lineage**. The
    run starts with ``chromosomes`` chromosomes of the given ``topology``, across which the
    ``initial_families`` founding genes are dealt **round-robin**; ``family_names=["toxin", …]`` additionally
    declares **named** families (remembered in ``result.family_names`` for ``result.has_family(node,
    "toxin")``), as in the family core; ``transfer_to`` / ``replacement`` / ``self_transfer`` behave
    as in the family core.

    The **chromosome tier** changes chromosome *number*: ``fission`` (split), ``fusion`` (merge — the
    reticulation), ``chromosome_origination`` (a de-novo replicon), ``chromosome_loss`` (a whole
    chromosome and its genes die; never the genome's last). Chromosomes carry identity — re-minted at
    every event that reshapes them — so ``chromosome_events`` is the true reticulating chromosome
    genealogy, rooted at the initial and de-novo originations. Deterministic given ``seed``.
    """
    tree = tree.complete_tree if isinstance(tree, SpeciesResult) else tree
    labels = _topologies(chromosomes, topology)
    n_initial_chrom = chromosomes
    dup = as_rate(duplication, default_scope=PerCopy)
    tra = as_rate(transfer, default_scope=PerCopy)
    los = as_rate(loss, default_scope=PerCopy)
    org = as_rate(origination, default_scope=PerLineage)
    inv = as_rate(inversion, default_scope=PerCopy)
    trp = as_rate(transposition, default_scope=PerCopy)
    trl = as_rate(translocation, default_scope=PerCopy)
    fis = as_rate(fission, default_scope=PerChromosome)
    fus = as_rate(fusion, default_scope=PerChromosome)
    cor = as_rate(chromosome_origination, default_scope=PerLineage)
    clo = as_rate(chromosome_loss, default_scope=PerChromosome)
    # this slice wires each event's default scope, OnTime (skyline), and ByFamily — the last with the
    # weight on the SEGMENT rather than on its starting gene (SPEC §6, and _pick_run_by_family). A
    # scope override or a clade/driven modifier is a later slice, so reject it rather than silently
    # mis-scale (see the family engine for the reasoning).
    for label, rate, want in (("duplication", dup, PerCopy), ("transfer", tra, PerCopy),
                              ("loss", los, PerCopy), ("origination", org, PerLineage),
                              ("inversion", inv, PerCopy), ("transposition", trp, PerCopy),
                              ("translocation", trl, PerCopy), ("fission", fis, PerChromosome),
                              ("fusion", fus, PerChromosome), ("chromosome_loss", clo, PerChromosome),
                              ("chromosome_origination", cor, PerLineage)):
        if not isinstance(rate.scope, want):
            raise ValueError(
                f"{label} has a {type(rate.scope).__name__} scope, but the ordered genome engine "
                f"wires only {want.__name__} for {label} this slice — scope overrides are a later slice."
            )
        for m in rate.modifiers:
            if isinstance(m, ByFamily) and label == "origination":
                raise ValueError(
                    "origination carries ByFamily, but origination is the rate at which families are "
                    "CREATED — when it is read there is no family yet to have drawn a factor for. "
                    "Put ByFamily on duplication, transfer, loss, inversion, transposition or "
                    "translocation, or use family_speed= for a family-wide tempo.")
            if isinstance(m, ByFamily) and not isinstance(rate.scope, PerCopy):
                raise ValueError(
                    f"{label} carries ByFamily on a {type(rate.scope).__name__} scope. A per-family "
                    f"weight has to reach the genes an event covers, so it is wired for the per-copy "
                    f"gene events only — not for the chromosome tier, which acts on whole replicons.")
            if not isinstance(m, (OnTime, ByFamily)):
                raise ValueError(
                    f"{label} carries {type(m).__name__}, which the ordered genome engine does not "
                    f"support yet — OnTime (skyline) and ByFamily are wired. Clade drift and driven "
                    f"rates are later slices."
                )
    # per-event extent distributions (segment size in genes); a bare number is the mean, None a single gene
    def _ext_spec(spec, label):
        """One event's extent (SPEC §6): ``base × modifiers``, no scope, in **genes** here. An extent
        takes the modifiers this resolution wires on a rate — ``OnTime`` — and they scale the size
        drawn. A driver is not among them: this engine wires no ``DrivenBy`` on its rates either, and
        for trait-driven rearrangement the nucleotide resolution is the one that has it."""
        e = as_extent(spec)
        for m in e.modifiers:
            if not isinstance(m, OnTime):
                raise ValueError(
                    f"{label} carries {type(m).__name__}, which the ordered genome engine does not "
                    f"wire on an extent — only OnTime (a skyline in time). For a trait-driven extent "
                    f"use --resolution nucleotide.")
        return e

    dup_ext, los_ext, tra_ext = (_ext_spec(duplication_extent, "duplication_extent"),
                                 _ext_spec(loss_extent, "loss_extent"),
                                 _ext_spec(transfer_extent, "transfer_extent"))
    inv_ext, trp_ext, trl_ext = (_ext_spec(inversion_extent, "inversion_extent"),
                                 _ext_spec(transposition_extent, "transposition_extent"),
                                 _ext_spec(translocation_extent, "translocation_extent"))
    if not 0.0 <= inversion_probability <= 1.0:
        raise ValueError(f"inversion_probability must be in [0, 1], got {inversion_probability!r}")
    if transfer_to == "distance":
        transfer_to = Distance()
    if transfer_to != "uniform" and not isinstance(transfer_to, Distance):
        raise ValueError(f"transfer_to must be 'uniform', 'distance', or Distance(decay=), got {transfer_to!r}")
    if isinstance(initial_families, bool) or not isinstance(initial_families, int) or initial_families < 0:
        raise ValueError(f"initial_families must be a non-negative integer, got {initial_families!r}")
    family_names = list(family_names) if family_names is not None else []
    for name in family_names:
        if not isinstance(name, str) or not name.strip():
            raise ValueError(f"family_names must be a list of non-empty family names (strings), got {name!r}")
    if len(set(family_names)) != len(family_names):
        raise ValueError(f"family names must be unique, got {family_names}")

    # The growth guard, as at the family resolution: duplication compounds, so a run whose rate sits
    # above its loss rate — or a family that drew a high ByFamily factor — multiplies without bound
    # unless something stops it. A segment may carry several families, and several copies of one, so
    # the run is refused when it would take *any* of them past the quota (see _run_over_cap).
    cap = resolve_max_family_size(max_family_size, len(tree.nodes))
    if family_speed is not None and not isinstance(family_speed, ByFamily):
        raise ValueError(
            f"family_speed must be a ByFamily(...) draw, got {type(family_speed).__name__}.")

    rng = np.random.default_rng(seed)
    copy_counter = 0
    family_counter = 0
    chrom_counter = 0

    def new_gene(family: int, strand: int) -> Gene:
        nonlocal copy_counter
        g = Gene(copy_counter, family, strand)
        copy_counter += 1
        return g

    # Per-family multipliers, drawn once when a family is minted and fixed for its whole life, exactly
    # as at the family resolution: family_speed scales every rate that family has (one draw), a
    # ByFamily on a single rate varies that rate alone (its own draw). What differs here is where the
    # weight lands — on the run an event covers, not on the gene it started from (SPEC §6). Empty
    # unless one of them is used, and a run without either is byte-identical to the plain path.
    fam_by = {"duplication": next((m for m in dup.modifiers if isinstance(m, ByFamily)), None),
              "transfer": next((m for m in tra.modifiers if isinstance(m, ByFamily)), None),
              "loss": next((m for m in los.modifiers if isinstance(m, ByFamily)), None),
              "inversion": next((m for m in inv.modifiers if isinstance(m, ByFamily)), None),
              "transposition": next((m for m in trp.modifiers if isinstance(m, ByFamily)), None),
              "translocation": next((m for m in trl.modifiers if isinstance(m, ByFamily)), None)}
    any_family = family_speed is not None or any(fam_by.values())
    fam_mult: dict[str, dict[int, float]] = {key: {} for key in fam_by}

    def new_family() -> int:
        nonlocal family_counter
        f = family_counter
        family_counter += 1
        if any_family:
            speed = family_speed.draw(rng) if family_speed is not None else 1.0
            for key, m in fam_by.items():
                fam_mult[key][f] = speed * (m.draw(rng) if m is not None else 1.0)
        return f

    def new_chromosome() -> int:
        nonlocal chrom_counter
        cid = chrom_counter
        chrom_counter += 1
        return cid

    depth = mean_root_to_tip(tree)  # timescale for Distance weighting (unused by "uniform")
    schedule = sorted((tree.nodes[i].end_time, i) for i in tree.nodes)  # (end_time, node_id)

    root = tree.nodes[tree.root]
    t = root.birth_time
    alive: list[int] = []
    gen: list[list[Chromosome]] = []
    pos: dict[int, int] = {}
    genomes: dict[int, tuple[Chromosome, ...]] = {}
    events: list[Event] = []
    event_positions: list[EventPosition] = []
    rearrangements: list[Inversion | Transposition | Translocation] = []
    chromosome_events: list[ChromosomeEvent] = []

    initial_chroms = []
    for label in labels:  # lay down the initial karyotype; each initial chromosome is a network root
        cid = new_chromosome()
        initial_chroms.append(Chromosome(cid, label, []))
        chromosome_events.append(ChromosomeEvent(t, "origination", root.id, (), (cid,)))
    # the crown's initial genome is logged like any other origination — each founding gene appended in turn —
    # so the position table is total over gene-content events and a replay of the root branch can
    # start from an empty karyotype (every other branch starts from its parent's gene_order rows)
    for i in range(initial_families):  # deal the founding genes round-robin across the chromosomes
        fam = new_family()
        chrom = initial_chroms[i % n_initial_chrom]
        chrom.genes.append(new_gene(fam, +1))
        events.append(Event(t, "origination", root.id, fam, chrom.genes[-1].id))
        event_positions.append(EventPosition(t, "origination", root.id, chrom.id,
                                             len(chrom.genes) - 1, 1, family=fam))
    named: dict[str, int] = {}  # a minted id per declared name, dealt round-robin after the anonymous ones
    for j, name in enumerate(family_names):
        fam = new_family()
        named[name] = fam
        chrom = initial_chroms[(initial_families + j) % n_initial_chrom]
        chrom.genes.append(new_gene(fam, +1))
        events.append(Event(t, "origination", root.id, fam, chrom.genes[-1].id))
        event_positions.append(EventPosition(t, "origination", root.id, chrom.id,
                                             len(chrom.genes) - 1, 1, family=fam))
    # the run's starting genome: a deep snapshot, so the live genome's events never reach it
    initial_genome = tuple(Chromosome(c.id, c.topology, list(c.genes)) for c in initial_chroms)
    enter(alive, gen, pos, root.id, initial_chroms)
    total_copies = initial_families + len(family_names)
    total_chromosomes = n_initial_chrom

    bar = progress_bar(len(schedule), "genomes", unit="branch", enabled=progress)
    si = 0
    while si < len(schedule):
        bar.to(si)
        n = total_copies
        k_alive = len(alive)
        ctx = {"copies": n, "lineages": k_alive, "chromosomes": total_chromosomes, "time": t}
        c = total_chromosomes
        can_xfer = n > 0 and (k_alive >= 2 or self_transfer)
        # A per-copy rate pools over genes, so with per-family weights the total is the unit rate
        # times those weights summed over the live genes — and the run must then be drawn with the
        # same weights, or the rate would say one thing and the picking another. Summed per lineage,
        # so the lineage pick can reuse them. On a circular chromosome ``Σ_s mean_w(s, m)`` is exactly
        # this sum for every run size, which is why no per-size term appears here (SPEC §6).
        fw = None
        if any_family:
            fw = {key: [sum(mult[g.family] for chrom in gen[k] for g in chrom.genes)
                        for k in range(k_alive)]
                  for key, mult in fam_mult.items()}
            one = {"copies": 1, "lineages": 1, "chromosomes": 1, "time": t}
            r_dup = dup.effective(**one) * sum(fw["duplication"]) if n else 0.0
            r_los = los.effective(**one) * sum(fw["loss"]) if n else 0.0
            r_tra = tra.effective(**one) * sum(fw["transfer"]) if can_xfer else 0.0
            r_inv = inv.effective(**one) * sum(fw["inversion"]) if n else 0.0
            r_trp = trp.effective(**one) * sum(fw["transposition"]) if n else 0.0
            r_trl = trl.effective(**one) * sum(fw["translocation"]) if n else 0.0
        else:
            r_dup = dup.effective(**ctx) if n else 0.0
            r_los = los.effective(**ctx) if n else 0.0
            r_tra = tra.effective(**ctx) if can_xfer else 0.0
            r_inv = inv.effective(**ctx) if n else 0.0                  # per copy (the run's start)
            r_trp = trp.effective(**ctx) if n else 0.0                  # per copy (the run's start)
            r_trl = trl.effective(**ctx) if n else 0.0                  # per copy; needs >=2 chromosomes
        r_org = org.effective(**ctx)                                    # per lineage
        r_fis = fis.effective(**ctx) if c else 0.0                      # per chromosome (the tier)
        r_fus = fus.effective(**ctx) if c else 0.0
        r_cor = cor.effective(**ctx)                                    # per lineage (de-novo replicon)
        r_clo = clo.effective(**ctx) if c else 0.0
        total = (r_dup + r_los + r_org + r_tra + r_inv + r_trp + r_trl
                 + r_fis + r_fus + r_cor + r_clo)

        next_species = schedule[si][0]
        horizon = min(next_species, dup.next_change(t), los.next_change(t), org.next_change(t),
                      tra.next_change(t), inv.next_change(t), trp.next_change(t), trl.next_change(t),
                      fis.next_change(t), fus.next_change(t), cor.next_change(t), clo.next_change(t))

        if total > 0.0:
            t_ev = t + float(rng.exponential(1.0 / total))
            if t_ev < horizon:  # a genome event fires before the alive set or a rate changes
                t = t_ev
                r = float(rng.random()) * total
                b_los = r_dup + r_los                    # cumulative bounds, in the firing order below
                b_org = b_los + r_org
                b_tra = b_org + r_tra
                b_inv = b_tra + r_inv
                b_trp = b_inv + r_trp
                b_trl = b_trp + r_trl
                b_fis = b_trl + r_fis
                b_fus = b_fis + r_fus
                b_cor = b_fus + r_cor                    # ... and the remainder (to total) is clo
                if r < r_dup:                            # every gene-level event acts on an extent
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "duplication", dup_ext, {"time": t})
                    if picked is not None:
                        k, ci, j, m = picked
                        if not _run_over_cap(gen[k], gen[k][ci], j, m, cap):
                            total_copies += _duplicate(gen[k][ci], j, m, tree.nodes[alive[k]], t,
                                                       events, event_positions, new_gene)
                elif r < b_los:
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "loss", los_ext, {"time": t})
                    if picked is not None:
                        k, ci, j, m = picked
                        total_copies -= _lose_at(gen[k][ci], j, m, tree.nodes[alive[k]], t, events,
                                                 event_positions)
                elif r < b_org:
                    k = int(rng.integers(k_alive))  # origination is per lineage: a uniform lineage
                    _originate(gen[k], tree.nodes[alive[k]], t, events, event_positions, new_gene,
                               new_family, rng)
                    total_copies += 1
                elif r < b_tra:
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "transfer", tra_ext, {"time": t})
                    if picked is not None:
                        kd, cdi, jd, m = picked
                        total_copies += _do_transfer(rng, tree, alive, gen, kd, cdi, jd, m, t, events,
                                                     event_positions, new_gene, transfer_to,
                                                     replacement, self_transfer, depth, cap)
                elif r < b_inv:
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "inversion", inv_ext, {"time": t})
                    if picked is not None:                # the run starts at a gene, so: per copy
                        k, ci, i0, m = picked
                        _invert(gen[k][ci], i0, m, tree.nodes[alive[k]], t, rearrangements)
                elif r < b_trp:
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "transposition", trp_ext, {"time": t})
                    if picked is not None:
                        k, ci, i0, m = picked
                        _transpose(gen[k][ci], i0, m, tree.nodes[alive[k]], t, rearrangements, rng,
                                   inversion_probability)
                elif r < b_trl:
                    picked = _pick_event_run(rng, gen, n, fw, fam_mult, "translocation", trl_ext, {"time": t})
                    if picked is not None:
                        k, ci, j, m = picked
                        _translocate(gen[k], ci, j, m, tree.nodes[alive[k]], t, rearrangements, rng,
                                     inversion_probability)
                elif r < b_fis:
                    k, ci = _pick_chromosome(rng, gen, c)
                    dc, dg = _fission(gen[k], ci, tree.nodes[alive[k]], t, chromosome_events,
                                      new_chromosome, rng)
                    total_chromosomes += dc
                    total_copies += dg
                elif r < b_fus:
                    k, ci = _pick_chromosome(rng, gen, c)
                    dc, dg = _fusion(gen[k], ci, tree.nodes[alive[k]], t, chromosome_events,
                                     new_chromosome, rng)
                    total_chromosomes += dc
                    total_copies += dg
                elif r < b_cor:
                    k = int(rng.integers(k_alive))  # chromosome origination is per lineage
                    dc, dg = _chromosome_originate(gen[k], tree.nodes[alive[k]], t, chromosome_events,
                                                   new_chromosome)
                    total_chromosomes += dc
                    total_copies += dg
                else:
                    k, ci = _pick_chromosome(rng, gen, c)
                    dc, dg = _chromosome_lose(gen[k], ci, tree.nodes[alive[k]], t, events,
                                              event_positions, chromosome_events)
                    total_chromosomes += dc
                    total_copies += dg
                continue

        if horizon == next_species:  # advance to the tree's next event(s); process the whole tie-batch
            t = next_species
            while si < len(schedule) and schedule[si][0] == t:
                i = schedule[si][1]
                g = gen[pos[i]]
                genomes[i] = tuple(Chromosome(c.id, c.topology, tuple(c.genes)) for c in g)  # freeze
                total_copies -= sum(len(c.genes) for c in g)
                total_chromosomes -= len(g)
                retire(alive, gen, pos, pos[i])
                node = tree.nodes[i]
                if node.children is not None:  # a speciation: re-mint every chromosome and gene id
                    child_genomes = {c: [] for c in node.children}
                    for pchrom in g:
                        dcids = []
                        for c in node.children:
                            dcid = new_chromosome()
                            dcids.append(dcid)
                            dgenes = []
                            for old in pchrom.genes:  # ZOMBI1: the gene ends and continues, fresh id
                                ng = new_gene(old.family, old.strand)
                                dgenes.append(ng)
                                events.append(Event(t, "speciation", c, old.family, ng.id, parent=old.id))
                            child_genomes[c].append(Chromosome(dcid, pchrom.topology, dgenes))
                        chromosome_events.append(
                            ChromosomeEvent(t, "speciation", node.id, (pchrom.id,), tuple(dcids)))
                    for c in node.children:
                        cg = child_genomes[c]
                        enter(alive, gen, pos, c, cg)
                        total_copies += sum(len(ch.genes) for ch in cg)
                        total_chromosomes += len(cg)
                si += 1
        else:
            t = horizon  # a skyline breakpoint: advance and re-evaluate the (now changed) rate

    bar.close()
    return OrderedGenomesResult(tree, genomes, events, rearrangements, chromosome_events, seed,
                                named, event_positions, initial_genome)

zombi2.genomes.simulate_genomes_nucleotide

simulate_genomes_nucleotide(tree, *, inversion=0.0, inversion_extent=50.0, translocation=0.0, translocation_extent=50.0, transposition=0.0, transposition_extent=50.0, inversion_probability=0.0, loss=0.0, loss_extent=50.0, duplication=0.0, duplication_extent=50.0, transfer=0.0, transfer_extent=50.0, transfer_to='uniform', self_transfer=False, origination=0.0, origination_extent=50.0, fission=0.0, fusion=0.0, chromosome_origination=0.0, chromosome_loss=0.0, chromosomes=1, root_length=1000, topology='circular', genes=0, gene_length=100, gff=None, fasta=None, trim_overlaps=False, seed=None, progress=False) -> NucleotideGenomesResult

Evolve a nucleotide genome along a species tree by inversion, translocation, transposition, loss, duplication, transfer, origination, and the number-changing chromosome tier. The run starts from a karyotype — chromosomes replicons, each its own source: an int N gives N equal replicons of root_length/topology, or pass a list of (length, topology) for heterogeneous sizes and shapes. Each lineage inherits a copy of its parent's karyotype at speciation, with every chromosome re-minted (the chromosome network), and evolves:

• inversion (per lineage) reverses a geometric-length (mean inversion_extent) arc of a length-weighted chromosome.
• translocation (per lineage) moves a geometric-length (mean translocation_extent) arc to a different chromosome; transposition (per lineage, mean transposition_extent) moves one within its chromosome. Both land inverted with probability inversion_probability, keep source coordinates, and are rearrangements, not edges.
• loss (per lineage) deletes a geometric-length (mean loss_extent) arc — an ancestry-changing event (a death), recorded in events. Never empties a chromosome.
• duplication (per lineage) copies a geometric-length (mean duplication_extent) arc in tandem — an ancestry-changing birth, recorded in events.
• transfer (per lineage) copies a geometric-length (mean transfer_extent) arc into a contemporaneous recipient (transfer_to: "uniform" or "distance" / a :class:Distance; self_transfer allows the donor itself) — a horizontal birth, additive (the donor keeps its copy). This is what needs the global timeline.
• origination (per lineage) lays down a new gene on a fresh source (geometric length, mean origination_extent) — a birth of a wholly new family, indivisible from birth.
• fission (per chromosome) splits a chromosome in two (a bifurcation); fusion (per chromosome) merges two chromosomes of the same topology (the reticulation). chromosome_origination (per lineage) adds a de-novo circular replicon (a plasmid, a network root) carrying one new gene; chromosome_loss (per chromosome) kills a whole chromosome (its material dies as a loss; never the last one) — a network leaf. All record a chromosome-network edge.

A chromosome never exists without a gene. A replicon is born with one, and any event that would strip a chromosome of its last gene — a loss, a translocation carrying it away, a fission splitting off a geneless half — simply does not happen. (Vacuous when no genes are declared.)

The engine runs a global-timeline Gillespie: all lineages alive at once evolve along one clock (every segmental event is per lineage — the rate says how often a lineage does it, the extent how much it touches, so a bigger genome does not get proportionally more events; the chromosome tier is per chromosome), so a transfer couples two contemporaries. With loss, the strong invariant weakens: every node carries a subset of the initial sequence (each ancestral position at most once, monotonically down every path); origination further adds fresh sources beyond the root. Deterministic given seed. (Transfer is always additive.)

Source code in zombi2/genomes/nucleotide.py

@without_cyclic_gc
def simulate_genomes_nucleotide(tree, *, inversion=0.0, inversion_extent=50.0, translocation=0.0,
                                translocation_extent=50.0, transposition=0.0, transposition_extent=50.0,
                                inversion_probability=0.0, loss=0.0, loss_extent=50.0, duplication=0.0,
                                duplication_extent=50.0, transfer=0.0, transfer_extent=50.0,
                                transfer_to="uniform", self_transfer=False, origination=0.0,
                                origination_extent=50.0, fission=0.0, fusion=0.0,
                                chromosome_origination=0.0, chromosome_loss=0.0, chromosomes=1,
                                root_length=1000, topology="circular", genes=0, gene_length=100,
                                gff=None, fasta=None, trim_overlaps=False, seed=None,
                                progress=False) -> NucleotideGenomesResult:
    """Evolve a nucleotide genome along a species tree by inversion, translocation, transposition,
    **loss**, **duplication**, **transfer**, **origination**, and the number-changing chromosome tier.
    The run starts from a **karyotype** — ``chromosomes`` replicons, each its own source: an int
    ``N`` gives ``N`` equal replicons of ``root_length``/``topology``, or pass a list of ``(length,
    topology)`` for heterogeneous **sizes and shapes**. Each lineage inherits a copy of its parent's
    karyotype at speciation, with **every chromosome re-minted** (the chromosome network), and evolves:

    - ``inversion`` (**per lineage**) reverses a geometric-length (mean ``inversion_extent``) arc
      of a length-weighted chromosome.
    - ``translocation`` (**per lineage**) moves a geometric-length (mean ``translocation_extent``)
      arc to a **different** chromosome; ``transposition`` (**per lineage**, mean
      ``transposition_extent``) moves one **within** its chromosome. Both land inverted with
      probability ``inversion_probability``, keep source coordinates, and are rearrangements, not edges.
    - ``loss`` (**per lineage**) deletes a geometric-length (mean ``loss_extent``) arc — an
      ancestry-**changing** event (a death), recorded in ``events``. Never empties a chromosome.
    - ``duplication`` (**per lineage**) copies a geometric-length (mean ``duplication_extent``) arc
      in tandem — an ancestry-**changing** *birth*, recorded in ``events``.
    - ``transfer`` (**per lineage**) copies a geometric-length (mean ``transfer_extent``) arc into a
      **contemporaneous recipient** (``transfer_to``: ``"uniform"`` or ``"distance"`` / a
      :class:`Distance`; ``self_transfer`` allows the donor itself) — a horizontal *birth*, additive
      (the donor keeps its copy). This is what needs the global timeline.
    - ``origination`` (**per lineage**) lays down a **new gene** on a fresh source (geometric length,
      mean ``origination_extent``) — a *birth* of a wholly new family, indivisible from birth.
    - ``fission`` (**per chromosome**) splits a chromosome in two (a **bifurcation**); ``fusion``
      (**per chromosome**) merges two chromosomes of the same topology (the **reticulation**).
      ``chromosome_origination`` (**per lineage**) adds a de-novo circular replicon (a plasmid, a
      network **root**) **carrying one new gene**; ``chromosome_loss`` (**per chromosome**) kills a
      whole chromosome (its material dies as a loss; never the last one) — a network **leaf**. All
      record a chromosome-network edge.

    **A chromosome never exists without a gene.** A replicon is born with one, and any event that
    would strip a chromosome of its last gene — a loss, a translocation carrying it away, a fission
    splitting off a geneless half — simply does not happen. (Vacuous when no genes are declared.)

    The engine runs a **global-timeline** Gillespie: all lineages alive at once evolve along one clock
    (every segmental event is **per lineage** — the rate says how often a lineage does it, the extent how
    much it touches, so a bigger genome does not get proportionally more events; the chromosome tier is
    per chromosome), so a transfer
    couples two contemporaries. With loss, the strong invariant weakens: every node carries a **subset**
    of the initial sequence (each ancestral position at most once, monotonically down every path);
    origination further adds fresh sources beyond the root. Deterministic given ``seed``. (Transfer is
    always additive.)"""
    tree = tree.complete_tree if isinstance(tree, SpeciesResult) else tree
    # Every rate takes the written form (SPEC §5). The scopes here are **per lineage** for the gene
    # events — the rate says how often a lineage does the event and the extent says how much DNA it
    # touches, so the number reads the same whatever the genome's size — and **per chromosome** for the
    # number-changing tier. A bare number therefore stays a bare number, and the scope is stated rather
    # than hardcoded, so it can be seen and (later) overridden.
    _scoped = (("inversion", inversion, PerLineage), ("translocation", translocation, PerLineage),
               ("transposition", transposition, PerLineage), ("loss", loss, PerLineage),
               ("duplication", duplication, PerLineage), ("transfer", transfer, PerLineage),
               ("origination", origination, PerLineage), ("fission", fission, PerChromosome),
               ("fusion", fusion, PerChromosome),
               ("chromosome_origination", chromosome_origination, PerLineage),
               ("chromosome_loss", chromosome_loss, PerChromosome))
    _rates: dict[str, Rate] = {}
    for label, spec, want in _scoped:
        if isinstance(spec, (int, float)) and not isinstance(spec, bool) and spec < 0:
            raise ValueError(f"{label} must be >= 0, got {spec}")
        r = as_rate(spec, default_scope=want)
        if not isinstance(r.scope, want):
            raise ValueError(
                f"{label} has a {type(r.scope).__name__} scope, but the nucleotide engine wires only "
                f"{want.__name__} for {label} this slice — scope overrides are a later slice.")
        for m in r.modifiers:
            if not isinstance(m, WIRED_MODIFIERS):
                raise ValueError(
                    f"{label} carries {type(m).__name__}, which the nucleotide genome engine does not "
                    f"support yet — only {', '.join(w.__name__ for w in WIRED_MODIFIERS)} is wired. "
                    f"Driving a nucleotide rate with a trait is the next slice.")
        _rates[label] = r
    def _as_bp_extent(spec, label):
        """An extent in base pairs (SPEC §6): ``base × modifiers``, no scope. A bare number *is* the
        mean, so ``500`` reads the same here as anywhere else.

        The base must be :class:`~zombi2.rates.distributions.Geometric` — this engine draws each arc's
        far end **directly from the genome's legal breakpoints** rather than drawing a size and
        clamping it, so an arbitrary shape would have to be re-weighted over that set instead of
        sampled. Refusing beats quietly approximating. The modifiers are the ones this resolution
        wires on a rate, and they scale the mean: an extent's modifier is read when an event fires, so
        it changes how much that event takes without touching any rate."""
        if isinstance(spec, (int, float)) and not isinstance(spec, bool) and spec < 1:
            raise ValueError(f"{label} must be >= 1 bp, got {spec}")
        e = as_extent(spec)
        if not isinstance(e.base, Geometric):
            raise ValueError(
                f"{label} has a {type(e.base).__name__} base, but the nucleotide engine wires a "
                f"geometric extent only — it draws each arc's far end directly from the legal "
                f"breakpoints, so another shape would have to be re-weighted over that set rather "
                f"than drawn. Pass a number (the mean in bp) or Geometric(mean=...).")
        for m in e.modifiers:
            if not isinstance(m, WIRED_MODIFIERS):
                raise ValueError(
                    f"{label} carries {type(m).__name__}, which the nucleotide genome engine does not "
                    f"support yet — an extent takes the same modifiers a rate does here "
                    f"({', '.join(w.__name__ for w in WIRED_MODIFIERS)}).")
        return e

    inversion_extent = _as_bp_extent(inversion_extent, "inversion_extent")
    translocation_extent = _as_bp_extent(translocation_extent, "translocation_extent")
    transposition_extent = _as_bp_extent(transposition_extent, "transposition_extent")
    loss_extent = _as_bp_extent(loss_extent, "loss_extent")
    duplication_extent = _as_bp_extent(duplication_extent, "duplication_extent")
    transfer_extent = _as_bp_extent(transfer_extent, "transfer_extent")
    origination_extent = _as_bp_extent(origination_extent, "origination_extent")
    _extents = {"inversion_extent": inversion_extent, "translocation_extent": translocation_extent,
                "transposition_extent": transposition_extent, "loss_extent": loss_extent,
                "duplication_extent": duplication_extent, "transfer_extent": transfer_extent,
                "origination_extent": origination_extent}
    if not 0.0 <= inversion_probability <= 1.0:
        raise ValueError(f"inversion_probability must be in [0, 1], got {inversion_probability}")
    if transfer_to == "distance":
        transfer_to = Distance()
    if transfer_to != "uniform" and not isinstance(transfer_to, Distance):
        raise ValueError(f"transfer_to must be 'uniform', 'distance', or Distance(decay=), "
                         f"got {transfer_to!r}")
    if isinstance(genes, bool) or not isinstance(genes, int) or genes < 0:
        raise ValueError(f"genes must be a non-negative integer, got {genes!r}")
    if genes and (isinstance(gene_length, bool) or not isinstance(gene_length, int) or gene_length < 1):
        raise ValueError(f"gene_length must be a positive integer, got {gene_length!r}")
    initial_sequence: dict[int, str] = {}               # {source: initial DNA}, empty unless a FASTA is given
    if gff is not None:                              # declared from a GFF: exact coordinates and names
        if genes:
            raise ValueError("pass either gff= or genes=, not both — a GFF already declares the genes")
        lengths, gff_genes = read_gff(gff, trim_overlaps=trim_overlaps)
        seqids = sorted(lengths)                     # a deterministic replicon order
        by_seqid: dict[str, list] = {sq: [] for sq in seqids}
        for gene in gff_genes:
            by_seqid[gene.seqid].append((gene.start, gene.end, gene.strand, gene.name))
        specs = [(_valid_length(lengths[sq]), topology) for sq in seqids]
        layouts = [by_seqid[sq] for sq in seqids]
        if fasta is not None:                        # the initial DNA, one record per replicon, by seqid
            seqs = read_fasta(fasta)
            if set(seqs) != set(seqids):
                raise ValueError(
                    f"the FASTA's records {sorted(seqs)} do not match the GFF's replicons "
                    f"{seqids} — every ##sequence-region needs exactly one > record, same id")
            for i, sq in enumerate(seqids):
                if len(seqs[sq]) != lengths[sq]:
                    raise ValueError(
                        f"replicon {sq!r} is {lengths[sq]} bp in the GFF but {len(seqs[sq])} bp in "
                        "the FASTA — the sequence must be exactly as long as its sequence-region")
                initial_sequence[i] = seqs[sq]
    else:
        if fasta is not None:
            raise ValueError("fasta= needs gff=: the FASTA's records are matched to the GFF's "
                             "replicons by id, so there is nothing to lay down without one")
        specs = _replicon_specs(chromosomes, root_length, topology)
        for _length, _top in specs:                  # the genes must fit; they need not leave a gap
            if genes and genes * gene_length > _length:
                raise ValueError(f"{genes} genes of {gene_length} bp do not fit in a {_length} bp "
                                 f"replicon")
        layouts = [_even_gene_intervals(length, genes, gene_length) for (length, _t) in specs]
    # Conditioning: a rate carrying DrivenBy reads a driver **per lineage**, so the rates stop being
    # one number for the whole live set and become one per lineage. Same machinery as the family
    # resolution — each source resolves once into a DriverTrajectory keyed by the shared species node
    # id, from a file or an in-memory trait result. With no driven rate this is empty and the loop
    # stays exactly the pooled one, so an uncoupled run is untouched.
    driven = {label: [m for m in r.modifiers if isinstance(m, DrivenBy)] for label, r in _rates.items()}
    ext_driven = {label: [m for m in e.modifiers if isinstance(m, DrivenBy)]
                  for label, e in _extents.items()}
    by_key: dict[object, object] = {}
    for mods in (*driven.values(), *ext_driven.values()):
        for m in mods:
            by_key.setdefault(m.key, m.source)
    resolved = {}
    if by_key:
        resolved = {key: resolve_driver(src, tree) for key, src in by_key.items()}
        # a mapping whose states never occur leaves every lineage on the default factor, so the run
        # would secretly be the uncoupled model — refuse it here, naming the driver
        for mods in (*driven.values(), *ext_driven.values()):
            for m in mods:
                label = m.source if isinstance(m.source, str) else f"<{type(m.source).__name__}>"
                check_mapping_fires(m.mapping, resolved[m.key].states(), source_label=label)
    # Only a driver on a **rate** makes the loop per-lineage and adds a Gillespie breakpoint. A driver
    # on an **extent** is read at the instant an event fires — it changes how much that event takes,
    # never how often one happens — so it deliberately stays out of `trajs`: no per-lineage rate
    # weights, no extra horizon steps. (SPEC §6.)
    _rate_keys = {m.key for mods in driven.values() for m in mods}
    trajs = {k: v for k, v in resolved.items() if k in _rate_keys}
    any_driven = bool(trajs)

    rates = _Rates(_rates["inversion"], _rates["translocation"], _rates["transposition"],
                   _rates["loss"], _rates["duplication"], _rates["transfer"], _rates["origination"],
                   _rates["fission"], _rates["fusion"], _rates["chromosome_origination"],
                   _rates["chromosome_loss"], inversion_extent,
                   translocation_extent, transposition_extent, loss_extent, duplication_extent,
                   transfer_extent, origination_extent, inversion_probability)
    depth = mean_root_to_tip(tree)                       # timescale for Distance weighting

    rng = np.random.default_rng(seed)
    chrom_counter = 0
    copy_counter = 0
    source_counter = len(specs)                          # de-novo sources continue past the initial sources

    def new_chrom_id() -> int:
        nonlocal chrom_counter
        cid = chrom_counter
        chrom_counter += 1
        return cid

    def new_copy() -> int:
        nonlocal copy_counter
        copy_counter += 1
        return copy_counter                             # copy ids start at 1 (0 = the unset sentinel)

    family_counter = 0

    def new_family() -> int:
        nonlocal family_counter
        family_counter += 1
        return family_counter                           # gene family ids start at 1 (0 = intergene)

    def new_source() -> int:
        nonlocal source_counter
        src = source_counter
        source_counter += 1
        return src

    genomes: dict[int, NucleotideGenome] = {}
    events: list[Origination | Loss | Duplication | Transfer | Speciation] = []
    rearrangements: list[Inversion | Translocation | Transposition] = []
    chromosome_events: list[ChromosomeEvent] = []
    root = tree.nodes[tree.root]
    schedule = sorted((tree.nodes[i].end_time, i) for i in tree.nodes)   # (end_time, node) in time order

    initial_chroms = []
    gene_spans: dict[int, tuple[int, int, int]] = {}
    gene_names: dict[str, int] = {}
    gene_strands: dict[int, int] = {}
    for source, ((length, top), intervals) in enumerate(zip(specs, layouts)):  # one source per replicon
        cid = new_chrom_id()
        cp = new_copy()                                 # ...and one initial copy lineage per replicon
        initial_chroms.append(Chromosome(cid, top, _initial_blocks(source, length, cp, intervals, new_family,
                                                             gene_spans, gene_names,
                                                             gene_strands)))
        chromosome_events.append(ChromosomeEvent(root.birth_time, "origination", root.id, (), (cid,)))
        events.append(Origination(root.birth_time, root.id, cid, cp, source, 0, length, initial=True))

    # the run's starting genome: a deep snapshot, so the live genome's events never reach it
    initial_genome = NucleotideGenome(
        [Chromosome(c.id, c.topology, [Block(b.source, b.start, b.end, b.strand, b.copy, b.gene)
                                       for b in c.blocks]) for c in initial_chroms])

    t = root.birth_time
    alive: list[int] = []                               # the live-lineage set (species._grow shape)
    gen: list[NucleotideGenome] = []
    pos: dict[int, int] = {}
    enter(alive, gen, pos, root.id, NucleotideGenome(initial_chroms))
    total_length = sum(c.length for c in initial_chroms)
    total_chromosomes = len(initial_chroms)

    bar = progress_bar(len(schedule), "genomes", unit="branch", enabled=progress)
    si = 0
    while si < len(schedule):
        bar.to(si)
        length, count, nlin = total_length, total_chromosomes, len(alive)
        can_xfer = nlin >= 2 or self_transfer
        # Each rate carries its own scope, so the count it is "per" comes from the context rather than
        # from a multiplication written here. The gene events are PER LINEAGE: the rate says how often
        # a lineage does the event and the extent says how much DNA it touches, so a bigger genome does
        # NOT get more events (that would double-count size and explode).
        ctx = {"copies": 0, "lineages": nlin, "chromosomes": count, "time": t}
        # A driven rate differs from lineage to lineage, so it is summed **over the living lineages**,
        # each read with its own driver value and its own chromosome count — and the weights are kept,
        # because the affected lineage must then be drawn with them too. An undriven rate stays pooled
        # (one .effective, uniform pick), so a run with no coupling is byte-identical to before.
        w: dict[str, list[float]] = {}
        if any_driven:
            drivers = [{key: trajs[key].value(alive[k], t) for key in trajs} for k in range(nlin)]
            for label, rate in _rates.items():
                if driven[label]:
                    w[label] = [rate.effective(copies=0, lineages=1,
                                               chromosomes=len(gen[k].chromosomes), time=t,
                                               drivers=drivers[k]) for k in range(nlin)]

        def _r(label, pooled, live=True):
            """The total for one event class: summed per-lineage when driven, pooled when not."""
            if not live:
                return 0.0
            return sum(w[label]) if label in w else pooled

        r_inv = _r("inversion", rates.inversion.effective(**ctx))
        r_trl = _r("translocation", rates.translocation.effective(**ctx))
        r_trp = _r("transposition", rates.transposition.effective(**ctx))
        r_los = _r("loss", rates.loss.effective(**ctx))
        r_dup = _r("duplication", rates.duplication.effective(**ctx))
        r_tra = _r("transfer", rates.transfer.effective(**ctx), live=can_xfer)
        r_org = _r("origination", rates.origination.effective(**ctx))
        r_fis = _r("fission", rates.fission.effective(**ctx))
        r_fus = _r("fusion", rates.fusion.effective(**ctx))
        r_cor = _r("chromosome_origination", rates.chromosome_origination.effective(**ctx))
        r_clo = _r("chromosome_loss", rates.chromosome_loss.effective(**ctx))
        total = (r_inv + r_trl + r_trp + r_los + r_dup + r_tra + r_org + r_fis + r_fus + r_cor + r_clo)
        next_species = schedule[si][0]
        # a skyline steps at a known time, so the race runs only to the next of those or the next
        # species event — whichever comes first — and the rates are re-read on the other side.
        horizon = min(next_species, rates.inversion.next_change(t), rates.translocation.next_change(t),
                      rates.transposition.next_change(t), rates.loss.next_change(t),
                      rates.duplication.next_change(t), rates.transfer.next_change(t),
                      rates.origination.next_change(t), rates.fission.next_change(t),
                      rates.fusion.next_change(t), rates.chromosome_origination.next_change(t),
                      rates.chromosome_loss.next_change(t))
        if any_driven:  # a driven rate also changes when its driver switches mid-branch — step there
            horizon = min(horizon, min((trajs[key].next_change(alive[k], t) for key in trajs
                                        for k in range(nlin)), default=math.inf))

        def _ext(label, k):
            """An extent's mean in bp for this event: the base mean, scaled by its modifiers read on
            the acting lineage. Undriven is the common case and costs nothing — the point of reading it
            here rather than in the rate loop is that an extent changes no rate, so it never had to be
            raced to."""
            e = _extents[label]
            if not e.is_driven:
                return e.base.mean
            return e.mean(time=t, drivers={key: resolved[key].value(alive[k], t) for key in resolved})

        def _pick(label, fallback=None):
            """The affected lineage: drawn by its own effective rate where that rate is driven — the
            same weights the total was summed with — and otherwise by the rate's own undriven rule,
            which is uniform for a per-lineage rate and by chromosome count for the tier."""
            ws = w.get(label)
            if ws:
                return weighted_index(rng, ws, sum(ws))
            return fallback() if fallback is not None else int(rng.integers(nlin))
        if total > 0.0:
            t_ev = t + float(rng.exponential(1.0 / total))
            if t_ev < horizon:                          # a genome event fires before the horizon
                t = t_ev
                r = float(rng.random()) * total
                b_trl = r_inv + r_trl
                b_trp = b_trl + r_trp
                b_los = b_trp + r_los
                b_dup = b_los + r_dup
                b_tra = b_dup + r_tra
                b_org = b_tra + r_org
                b_fis = b_org + r_fis
                b_fus = b_fis + r_fus
                b_cor = b_fus + r_cor
                if r < r_inv:
                    k = _pick("inversion")
                    _do_inversion(gen[k], alive[k], t, _ext("inversion_extent", k), rng, rearrangements)
                elif r < b_trl:
                    k = _pick("translocation")
                    _do_translocation(gen[k], alive[k], t, _ext("translocation_extent", k),
                                      rates.inversion_probability, rng, rearrangements)
                elif r < b_trp:
                    k = _pick("transposition")
                    _do_transposition(gen[k], alive[k], t, _ext("transposition_extent", k),
                                      rates.inversion_probability, rng, rearrangements)
                elif r < b_los:
                    k = _pick("loss")
                    total_length += _do_loss(gen[k], alive[k], t, _ext("loss_extent", k), rng, events)
                elif r < b_dup:
                    k = _pick("duplication")
                    total_length += _do_duplication(gen[k], alive[k], t, _ext("duplication_extent", k),
                                                    rng, events, new_copy)
                elif r < b_tra:
                    kd = _pick("transfer")
                    total_length += _do_transfer(rng, tree, alive, gen, kd, t, _ext("transfer_extent", kd),
                                                 transfer_to, self_transfer, depth, events, new_copy)
                elif r < b_org:
                    k = _pick("origination")            # per lineage; weighted when driven
                    total_length += _do_origination(gen[k], alive[k], t, _ext("origination_extent", k),
                                                    rng, events, new_source, new_copy,
                                                    new_family, gene_spans, gene_strands)
                elif r < b_fis:
                    k = _pick("fission", lambda: _pick_lineage_by_chromosomes(rng, gen, count))
                    total_chromosomes += _do_fission(gen[k], alive[k], t, rng, chromosome_events,
                                                     new_chrom_id)
                elif r < b_fus:
                    k = _pick("fusion", lambda: _pick_lineage_by_chromosomes(rng, gen, count))
                    total_chromosomes += _do_fusion(gen[k], alive[k], t, rng, chromosome_events,
                                                    new_chrom_id)
                elif r < b_cor:
                    k = _pick("chromosome_origination")  # per lineage; weighted when driven
                    dc, dl = _do_chromosome_origination(
                        gen[k], alive[k], t, _ext("origination_extent", k), rng, events,
                        chromosome_events, new_chrom_id, new_source, new_copy, new_family,
                        gene_spans, gene_strands)
                    total_chromosomes += dc
                    total_length += dl
                else:
                    k = _pick("chromosome_loss", lambda: _pick_lineage_by_chromosomes(rng, gen, count))
                    dc, dl = _do_chromosome_loss(gen[k], alive[k], t, rng, events, chromosome_events)
                    total_chromosomes += dc
                    total_length += dl
                continue

        t = horizon
        if horizon < next_species:                      # a rate stepped, not a species event: re-read
            continue                                    # the rates on the other side and race again
        while si < len(schedule) and schedule[si][0] == t:   # process the whole tie-batch
            i = schedule[si][1]
            g = gen[pos[i]]
            genomes[i] = g                              # freeze: the lineage retires, never mutated again
            total_length -= g.length
            total_chromosomes -= len(g.chromosomes)
            retire(alive, gen, pos, pos[i])
            node = tree.nodes[i]
            if node.children is not None:              # a speciation: re-mint into the daughters
                for c, cg in _speciate(node, g, new_chrom_id, new_copy, events,
                                       chromosome_events).items():
                    enter(alive, gen, pos, c, cg)
                    total_length += cg.length
                    total_chromosomes += len(cg.chromosomes)
            si += 1
    bar.close()
    return NucleotideGenomesResult(tree, genomes, events, rearrangements, chromosome_events, seed,
                                  gene_spans, gene_names, gene_strands, initial_genome, initial_sequence)

zombi2.genomes.FamilyGenomesResult dataclass

FamilyGenomesResult(complete_tree: Tree, genomes: dict[int, tuple[GeneCopy, ...]], events: list[Event], seed: int | None, family_names: dict[str, int] = dict(), initial_genome: tuple[GeneCopy, ...] = ())

What simulate_genomes_family returns: the complete_tree it ran on, the final genomes at every node (extant and extinct), the events log (the compact source of truth), and the seed. The observed genomes are the extant tips — {n.id: genomes[n.id] for n in complete_tree.extant()}. The phyletic profiles are derived from those tips on access, and write materialises the chosen outputs to disk.

profiles cached property

profiles: Profiles

The phyletic profiles — each gene family's copy count in each extant species — derived from the observed genomes (the classic comparative-genomics matrix). See :mod:.profiles.

gene_trees cached property

gene_trees: dict[int, GeneTree]

{family id: GeneTree} — each family's true genealogy inside the complete tree, derived from the event log. Each GeneTree exposes .complete and .extant. See :mod:.gene_trees.

family_counts

family_counts(node_id: int) -> collections.Counter

A multiset view of one node's genome: family id → copy count.

Source code in zombi2/genomes/family.py

def family_counts(self, node_id: int) -> collections.Counter:
    """A multiset view of one node's genome: ``family id → copy count``."""
    return collections.Counter(c.family for c in self.genomes[node_id])

has_family

has_family(node_id: int, name: str) -> bool

Whether the named family name (declared via family_names=) is present — has ≥ 1 copy — in the genome at node_id. The presence signal a joint DrivenBy("genomes:<name>", …) reads.

Source code in zombi2/genomes/family.py

def has_family(self, node_id: int, name: str) -> bool:
    """Whether the named family ``name`` (declared via ``family_names=``) is present — has ≥ 1 copy — in
    the genome at ``node_id``. The presence signal a joint ``DrivenBy("genomes:<name>", …)`` reads."""
    if name not in self.family_names:
        raise KeyError(f"no named family {name!r}; declared families are {sorted(self.family_names)}")
    fid = self.family_names[name]
    return any(c.family == fid for c in self.genomes[node_id])

write

write(directory, outputs=('events', 'profiles', 'genomes', 'initial_genome', 'gene_trees')) -> None

Materialise chosen outputs to directory (created if needed):

• "events" → genome_events.tsv, the event log (the source of truth).
• "profiles" → profiles.tsv, the family × extant-species copy-count matrix.
• "genomes" → genomes.tsv, every node's gene content, one row per gene copy — ancestors included, where profiles.tsv counts only the extant tips.
• "initial_genome" → initial_genome.tsv, the genome the run started with. Its own file, not a row in genomes.tsv, because it belongs to no node: it sits at the start of the root branch, and every lineage in that table is a node at the end of one.
• "gene_trees" → gene_tree_fam<family>_{complete,extant}.nwk, each family's true genealogy. A family with no surviving copy writes no _extant file.

Source code in zombi2/genomes/family.py

def write(self, directory,
          outputs=("events", "profiles", "genomes", "initial_genome", "gene_trees")) -> None:
    """Materialise chosen ``outputs`` to ``directory`` (created if needed):

    - ``"events"`` → ``genome_events.tsv``, the event log (the source of truth).
    - ``"profiles"`` → ``profiles.tsv``, the family × extant-species copy-count matrix.
    - ``"genomes"`` → ``genomes.tsv``, every node's gene content, one row per gene copy —
      **ancestors included**, where ``profiles.tsv`` counts only the extant tips.
    - ``"initial_genome"`` → ``initial_genome.tsv``, the genome the run started with. Its own
      file, not a row in ``genomes.tsv``, because it belongs to no node: it sits at the start of
      the root branch, and every ``lineage`` in that table is a node at the end of one.
    - ``"gene_trees"`` → ``gene_tree_fam<family>_{complete,extant}.nwk``, each family's true
      genealogy. A family with no surviving copy writes no ``_extant`` file.
    """
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    if "events" in outputs:
        (d / "genome_events.tsv").write_text(events_tsv(self.events))
    if "profiles" in outputs:
        (d / "profiles.tsv").write_text(self.profiles.to_tsv())
    if "genomes" in outputs:
        (d / "genomes.tsv").write_text(self._genomes_tsv())
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv())
    if "gene_trees" in outputs:
        write_gene_trees(self.gene_trees, d)

zombi2.genomes.OrderedGenomesResult dataclass

OrderedGenomesResult(complete_tree: Tree, genomes: dict[int, tuple[Chromosome, ...]], events: list[Event], rearrangements: list[Inversion | Transposition | Translocation], chromosome_events: list[ChromosomeEvent], seed: int | None, family_names: dict[str, int] = dict(), event_positions: list[EventPosition] = list(), initial_genome: tuple[Chromosome, ...] = ())

What :func:simulate_genomes_ordered returns: the complete_tree it ran on, the final genomes at every node as tuples of :class:Chromosome\ s, the shared gene-genealogy events log, the rearrangements (inversions) and chromosome_events (the chromosome genealogy) logs, and the seed. The observed genomes are the extant tips; profiles and gene_trees are derived from the (position-blind) genealogy exactly as for the family core; gene_order reads a node's layout, and write materialises the chosen outputs.

profiles cached property

profiles: Profiles

The phyletic profiles — each gene family's copy count in each extant species — derived from the observed genomes, flattening across chromosomes (position does not enter). See :mod:.profiles.

gene_trees cached property

gene_trees: dict[int, GeneTree]

{family id: GeneTree} — each family's true genealogy inside the complete tree, derived from the (position-blind) event log exactly as for the family core. See :mod:.gene_trees.

family_counts

family_counts(node_id: int) -> collections.Counter

A multiset view of one node's genome: family id → copy count (across all chromosomes).

Source code in zombi2/genomes/ordered.py

def family_counts(self, node_id: int) -> collections.Counter:
    """A multiset view of one node's genome: ``family id → copy count`` (across all chromosomes)."""
    return collections.Counter(g.family for chrom in self.genomes[node_id] for g in chrom.genes)

has_family

has_family(node_id: int, name: str) -> bool

Whether the named family name (declared via family_names=) has ≥ 1 copy in the genome at node_id (across all chromosomes).

Source code in zombi2/genomes/ordered.py

def has_family(self, node_id: int, name: str) -> bool:
    """Whether the named family ``name`` (declared via ``family_names=``) has ≥ 1 copy in the genome at
    ``node_id`` (across all chromosomes)."""
    if name not in self.family_names:
        raise KeyError(f"no named family {name!r}; declared families are {sorted(self.family_names)}")
    fid = self.family_names[name]
    return any(g.family == fid for chrom in self.genomes[node_id] for g in chrom.genes)

gene_order

gene_order(node_id: int) -> list[tuple[int, int, int, int, int]]

One node's layout as (chromosome, position, strand, family, gene id) rows, chromosome by chromosome and left to right within each — the ordered analogue of family_counts.

Source code in zombi2/genomes/ordered.py

def gene_order(self, node_id: int) -> list[tuple[int, int, int, int, int]]:
    """One node's layout as ``(chromosome, position, strand, family, gene id)`` rows, chromosome
    by chromosome and left to right within each — the ordered analogue of ``family_counts``."""
    return [(chrom.id, pos, g.strand, g.family, g.id)
            for chrom in self.genomes[node_id] for pos, g in enumerate(chrom.genes)]

write

write(directory, outputs=('events', 'profiles', 'gene_order', 'initial_genome', 'gene_trees', 'chromosome_events')) -> None

Materialise chosen outputs to directory (created if needed):

• "events" → genome_events.tsv, the run's whole history in one time-ordered table: the gene genealogy, where each event happened, and the ancestry-neutral rearrangements. With gene_order this is enough to replay the run.
• "profiles" → profiles.tsv, the family × extant-species copy-count matrix.
• "gene_order" → gene_order.tsv, every node's layout (one row per gene), ancestors included — so a branch's rearrangements can be replayed from its parent's genome.
• "initial_genome" → initial_genome.tsv, the layout the run started with. Its own file, not a row in gene_order.tsv, because it belongs to no node: it sits at the start of the root branch, and every lineage in that table is a node at the end of one.
• "chromosome_events" → chromosome_events.tsv, the chromosome genealogy edges. The one log kept apart: it is a network over chromosome ids, with list-valued parents and children, joined on a different key from everything above.
• "gene_trees" → gene_tree_fam<family>_{complete,extant}.nwk, each family's true genealogy — unchanged from the family resolution, position being orthogonal to it.

Source code in zombi2/genomes/ordered.py

def write(self, directory,
          outputs=("events", "profiles", "gene_order", "initial_genome", "gene_trees",
                   "chromosome_events")) -> None:
    """Materialise chosen ``outputs`` to ``directory`` (created if needed):

    - ``"events"`` → ``genome_events.tsv``, the run's whole history in one time-ordered table:
      the gene genealogy, **where** each event happened, and the ancestry-neutral
      rearrangements. With ``gene_order`` this is enough to replay the run.
    - ``"profiles"`` → ``profiles.tsv``, the family × extant-species copy-count matrix.
    - ``"gene_order"`` → ``gene_order.tsv``, every node's layout (one row per gene), ancestors
      included — so a branch's rearrangements can be replayed from its parent's genome.
    - ``"initial_genome"`` → ``initial_genome.tsv``, the layout the run started with. Its own
      file, not a row in ``gene_order.tsv``, because it belongs to no node: it sits at the start
      of the root branch, and every ``lineage`` in that table is a node at the end of one.
    - ``"chromosome_events"`` → ``chromosome_events.tsv``, the chromosome genealogy edges. The
      one log kept apart: it is a network over chromosome **ids**, with list-valued parents and
      children, joined on a different key from everything above.
    - ``"gene_trees"`` → ``gene_tree_fam<family>_{complete,extant}.nwk``, each family's true
      genealogy — unchanged from the family resolution, position being orthogonal to it.
    """
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    if "events" in outputs:
        (d / "genome_events.tsv").write_text(
            _events_tsv(self.events, self.event_positions, self.rearrangements))
    if "profiles" in outputs:
        (d / "profiles.tsv").write_text(self.profiles.to_tsv())
    if "gene_order" in outputs:
        (d / "gene_order.tsv").write_text(self._gene_order_tsv())
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv())
    if "chromosome_events" in outputs:
        (d / "chromosome_events.tsv").write_text(chromosome_events_tsv(self.chromosome_events))
    if "gene_trees" in outputs:
        write_gene_trees(self.gene_trees, d)

zombi2.genomes.NucleotideGenomesResult dataclass

NucleotideGenomesResult(complete_tree: Tree, genomes: dict[int, NucleotideGenome], events: list[Origination | Loss | Duplication | Transfer | Speciation], rearrangements: list[Inversion | Translocation | Transposition], chromosome_events: list[ChromosomeEvent], seed: int | None, gene_spans: dict[int, tuple[int, int, int]] = dict(), gene_names: dict[str, int] = dict(), gene_strands: dict[int, int] = dict(), initial_genome: NucleotideGenome = (lambda: NucleotideGenome([]))(), initial_sequence: dict[int, str] = dict())

What :func:simulate_genomes_nucleotide returns: the complete_tree it ran on, the final nucleotide genomes (karyotypes) at every node, the copy-lineage genealogy events (origination, loss, duplication, transfer, speciation — carrying the copy ids the gene-tree recovery reads), the ancestry-neutral rearrangements (inversion, translocation), the chromosome_events (the chromosome network), and the seed. mosaic / trace_back / ancestry read a node's genome.

root_blocks property

root_blocks: list[tuple[int, int, int]]

The recovered root partition: (source, start, end) for each maximal never-cut interval that some node still carries — one per :attr:block_trees entry (by index).

Cut at the breakpoints of every node's genome, not only the extant leaves', which is what lets any node be reconstructed (:meth:assembly) rather than the survivors alone.

block_trees property

block_trees: dict[int, GeneTree]

{root-block index: GeneTree} — a tree for every recovered root block, spacer as well as genes, keyed by its index in :attr:root_blocks.

:attr:gene_trees covers the declared genes; this covers the whole genome. A block never splits, so its size is fixed and its genealogy is in the event log just as a gene's is — the recovery is the same one, pointed at every block instead of a chosen few. That is what makes an ancestral genome reconstructable at any node rather than only at the loci you declared.

A gene's tree here has the same topology and branch lengths as its :attr:gene_trees one, but not the same g<id> leaf labels: segment ids are handed out as the recovery walks its targets, and walking every block numbers them differently from walking three. Use one accessor or the other within a piece of analysis — they are the same genealogy under different names.

gene_trees property

gene_trees: dict[int, GeneTree]

{family: GeneTree} — the recovered gene trees.

With genes declared, one tree per gene, keyed by its gene family id (see :attr:gene_spans); the intergenic root-blocks keep their block ancestry in the log but are not built into trees. With no genes declared the whole genome is one big intergene, so every recovered root-block is a family in its own right and the key is its index in :attr:root_blocks. Every node votes on the partition, so a gene surviving only in lineages that died still gets a tree — a complete one, with no extant tree to go with it. Only a gene lost from every node has no root-block and no tree.

block_of

block_of(family: int) -> int

The index in :attr:root_blocks of the block a declared gene family occupies — the join between the two numbering schemes this resolution has.

They are both plain ints over overlapping ranges, so mixing them up is silent: gene_spans and gene_trees are keyed by gene family id, while root_blocks, block_trees and everything a sequence run produces here are keyed by block index (every block evolves, and spacer has no family). block_of is how you get from a gene to its sequences::

r.alignments[g.block_of(g.gene_names["dnaA"])]     # that gene's alignment

Raises KeyError for a family that was never declared, and LookupError for one declared but surviving nowhere at all — it has no recovered block, so there is nothing to point at. The reverse lookup is one line: {span: fam for fam, span in g.gene_spans.items()} read at root_blocks[i].

Source code in zombi2/genomes/nucleotide.py

def block_of(self, family: int) -> int:
    """The index in :attr:`root_blocks` of the block a declared **gene family** occupies — the join
    between the two numbering schemes this resolution has.

    They are both plain ints over overlapping ranges, so mixing them up is silent: ``gene_spans``
    and ``gene_trees`` are keyed by **gene family id**, while ``root_blocks``, ``block_trees`` and
    everything a sequence run produces here are keyed by **block index** (every block evolves, and
    spacer has no family). ``block_of`` is how you get from a gene to its sequences::

        r.alignments[g.block_of(g.gene_names["dnaA"])]     # that gene's alignment

    Raises ``KeyError`` for a family that was never declared, and ``LookupError`` for one declared
    but surviving nowhere at all — it has no recovered block, so there is nothing to point at.
    The reverse lookup is one line: ``{span: fam for fam, span in g.gene_spans.items()}`` read at
    ``root_blocks[i]``."""
    span = self.gene_spans[family]                    # KeyError: never declared
    try:
        return self.root_blocks.index(span)
    except ValueError:
        raise LookupError(
            f"gene family {family} spans {span} but has no recovered root block — no node in the "
            "tree still carries it, so nothing was reconstructed for it. gene_trees leaves such a "
            "family out for the same reason.") from None

assembly

assembly(node_id: int) -> dict[int, list[tuple[int, int, int]]]

How this node's genome is built out of the recovered root blocks: {chromosome id: [(block, gene, strand), …]} in physical order, where block indexes :attr:root_blocks, gene is the gene id that block's tree gives this node's copy (the g<id> label in :attr:block_trees), and strand is +1 read forward or -1 reverse-complemented.

To reconstruct a genome: pair each piece with its block's evolved sequence, flip the -1\ s, and concatenate. The sequence level does exactly that; nothing here knows about letters. Every node works — an extinct leaf and the root as readily as a surviving tip — which is what makes the whole history recoverable rather than only its leaves. :meth:initial_assembly does the same for the genome the run started with.

A piece is always a whole block, never part of one, because every node votes on where the partition is cut (see :func:_root_block_partition): this node's own breakpoints are all in it, so each of its blocks is a whole number of root blocks. What a block is cut into is one piece per root block it spans — and on a reversed block those come out in descending coordinate order, since physical order runs down the source.

Source code in zombi2/genomes/nucleotide.py

def assembly(self, node_id: int) -> dict[int, list[tuple[int, int, int]]]:
    """How this node's genome is built out of the recovered root blocks:
    ``{chromosome id: [(block, gene, strand), …]}`` in **physical order**, where ``block`` indexes
    :attr:`root_blocks`, ``gene`` is the gene id that block's tree gives this node's copy (the
    ``g<id>`` label in :attr:`block_trees`), and ``strand`` is ``+1`` read forward or ``-1``
    reverse-complemented.

    To reconstruct a genome: pair each piece with its block's evolved sequence, flip the
    ``-1``\\ s, and concatenate. The sequence level does exactly that; nothing here knows about
    letters. **Every** node works — an extinct leaf and the root as readily as a surviving tip —
    which is what makes the whole history recoverable rather than only its leaves.
    :meth:`initial_assembly` does the same for the genome the run started with.

    A piece is always a **whole** block, never part of one, because every node votes on where the
    partition is cut (see :func:`_root_block_partition`): this node's own breakpoints are all in
    it, so each of its blocks is a whole number of root blocks. What a block *is* cut into is one
    piece per root block it spans — and on a reversed block those come out in descending
    coordinate order, since physical order runs *down* the source."""
    tips = self._recover_blocks()[2]
    blocks = self.root_blocks
    what = node_label(node_id)
    out = {}
    for cid, pieces in self._pieces(self.genomes[node_id], what).items():
        named = []
        for (i, copy, strand) in pieces:
            gene = tips.get((i, copy), _MISSING)
            if gene is _MISSING or gene is None:
                raise AssertionError(                            # a guard — see the class docstring
                    f"{what} carries {blocks[i]} under copy lineage {copy}, but that block's "
                    + ("genealogy has no such copy" if gene is _MISSING
                       else "genealogy ends that copy in a loss")
                    + " — the event log and the genomes disagree")
            named.append((i, gene, strand))
        out[cid] = named
    return out

initial_assembly

initial_assembly() -> dict[int, list[tuple[int, int]]]

:meth:assembly for :attr:initial_genome: {chromosome id: [(block, strand), …]}.

No gene id here, unlike :meth:assembly, and that is the honest shape rather than a saving. The initial genome sits at the start of the root branch, before any event, so each of its blocks has exactly one sequence — the founding draw the sequence level records as founding[block] — and there is no copy to disambiguate. A gene id would in fact be wrong here: the one :meth:assembly gives is the last gene a copy held, and for an initial copy that is at the far end of the stem. A loss on the stem can even end it, which is the same thing said louder.

Source code in zombi2/genomes/nucleotide.py

def initial_assembly(self) -> dict[int, list[tuple[int, int]]]:
    """:meth:`assembly` for :attr:`initial_genome`: ``{chromosome id: [(block, strand), …]}``.

    No gene id here, unlike :meth:`assembly`, and that is the honest shape rather than a saving.
    The initial genome sits at the **start** of the root branch, before any event, so each of its
    blocks has exactly one sequence — the founding draw the sequence level records as
    ``founding[block]`` — and there is no copy to disambiguate. A gene id would in fact be *wrong*
    here: the one :meth:`assembly` gives is the **last** gene a copy held, and for an initial copy
    that is at the far end of the stem. A loss on the stem can even end it, which is the same
    thing said louder."""
    return {cid: [(i, strand) for (i, _copy, strand) in pieces]
            for cid, pieces in self._pieces(self.initial_genome, "the initial genome").items()}

write

write(directory, outputs=('events', 'genes', 'blocks', 'initial_genome', 'initial_sequence', 'gene_trees', 'chromosome_events', 'gff', 'bed')) -> None

Materialise chosen outputs to directory (created if needed):

• "events" → genome_events.tsv, the run's whole history in one time-ordered table: the copy-lineage genealogy and the ancestry-neutral rearrangements. One row per ancestral interval an event touched, so an event that spanned several blocks writes several rows sharing a time and kind.
• "blocks" → blocks.tsv, every node's genome as its block mosaic (ancestors included, as for the ordered resolution's gene_order). The one big file here: blocks are not kept maximal during a run, so a rearrangement-heavy genome carries far more of them than it has distinct ancestral runs, and this grows with their number × every node.
• "genes" → genes.tsv, the declared genes and where they sit in root coordinates. Header-only for a run that declared none.
• "initial_genome" → initial_genome.tsv, the block mosaic the run started with. Its own file, not a row in blocks.tsv, because it belongs to no node: it sits at the start of the root branch, and every lineage in that table is a node at the end of one.
• "chromosome_events" → chromosome_events.tsv, the chromosome network's edges. The one log kept apart: it is a network over chromosome ids, with list-valued parents and children, joined on a different key from everything above.
• "gene_trees" → gene_tree_fam<family>_{complete,extant}.nwk, one recovered genealogy per family some node still carries; the _extant file only where the family has a surviving copy.
• "initial_sequence" → initial_sequence.fasta, the initial DNA the run was given (fasta=), one >source<n> record per replicon. Written only when a FASTA was supplied — it is what lets a separate zombi2 sequences run found its blocks from the real sequence.
• "gff" → genome_<lineage>.gff, that genome's genes, in its own coordinates: the annotation to read beside the sequence level's genome_<lineage>.fasta.
• "bed" → genome_<lineage>.bed, that genome's blocks — every piece, spacer included, named by the ancestral interval it descends from. The ancestry as a browser track.

Both name their sequences <lineage>_chr<c>, exactly as the FASTA records are named, so a genome and its annotation join without renaming anything. Written for every node and for the initial genome, so there are two files per genome.

Source code in zombi2/genomes/nucleotide.py

def write(self, directory,
          outputs=("events", "genes", "blocks", "initial_genome", "initial_sequence",
                   "gene_trees", "chromosome_events", "gff", "bed")) -> None:
    """Materialise chosen ``outputs`` to ``directory`` (created if needed):

    - ``"events"`` → ``genome_events.tsv``, the run's whole history in one time-ordered table:
      the copy-lineage genealogy and the ancestry-neutral rearrangements. One row per
      **ancestral interval** an event touched, so an event that spanned several blocks writes
      several rows sharing a ``time`` and ``kind``.
    - ``"blocks"`` → ``blocks.tsv``, every node's genome as its block mosaic (ancestors
      included, as for the ordered resolution's ``gene_order``). The one big file here: blocks
      are not kept maximal during a run, so a rearrangement-heavy genome carries far more of
      them than it has distinct ancestral runs, and this grows with their number × every node.
    - ``"genes"`` → ``genes.tsv``, the declared genes and where they sit in root coordinates.
      Header-only for a run that declared none.
    - ``"initial_genome"`` → ``initial_genome.tsv``, the block mosaic the run started with. Its
      own file, not a row in ``blocks.tsv``, because it belongs to no node: it sits at the start
      of the root branch, and every ``lineage`` in that table is a node at the end of one.
    - ``"chromosome_events"`` → ``chromosome_events.tsv``, the chromosome network's edges. The
      one log kept apart: it is a network over chromosome **ids**, with list-valued parents and
      children, joined on a different key from everything above.
    - ``"gene_trees"`` → ``gene_tree_fam<family>_{complete,extant}.nwk``, one recovered
      genealogy per family some node still carries; the ``_extant`` file only where the family
      has a surviving copy.
    - ``"initial_sequence"`` → ``initial_sequence.fasta``, the initial DNA the run was given (``fasta=``),
      one ``>source<n>`` record per replicon. Written only when a FASTA was supplied — it is what
      lets a separate ``zombi2 sequences`` run found its blocks from the real sequence.
    - ``"gff"`` → ``genome_<lineage>.gff``, that genome's **genes**, in its own coordinates: the
      annotation to read beside the sequence level's ``genome_<lineage>.fasta``.
    - ``"bed"`` → ``genome_<lineage>.bed``, that genome's **blocks** — every piece, spacer
      included, named by the ancestral interval it descends from. The ancestry as a browser track.

    Both name their sequences ``<lineage>_chr<c>``, exactly as the FASTA records are named, so a
    genome and its annotation join without renaming anything. Written for every node and for the
    initial genome, so there are two files per genome.
    """
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    if "events" in outputs:
        (d / "genome_events.tsv").write_text(
            _nucleotide_events_tsv(self.events, self.rearrangements))
    if "blocks" in outputs:
        (d / "blocks.tsv").write_text(self._blocks_tsv())
    if "genes" in outputs:
        (d / "genes.tsv").write_text(self._genes_tsv())
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv())
    if "chromosome_events" in outputs:
        (d / "chromosome_events.tsv").write_text(chromosome_events_tsv(self.chromosome_events))
    if "gene_trees" in outputs:
        write_gene_trees(self.gene_trees, d)
    if "initial_sequence" in outputs and self.initial_sequence:
        (d / "initial_sequence.fasta").write_text(
            "".join(f">source{src}\n{self.initial_sequence[src]}\n"
                    for src in sorted(self.initial_sequence)))
    for token, ext, render in (("gff", "gff", self._gff), ("bed", "bed", self._bed)):
        if token in outputs:
            for label, genome in self._every_genome():
                (d / f"genome_{label}.{ext}").write_text(render(label, genome))

zombi2.genomes.GeneTree dataclass

GeneTree(family: int, complete: GeneNode, origination: float)

One gene family's true genealogy. complete is the whole tree (lost and extinct-species lineages included); extant is it pruned to the genes surviving at the extant tips (degree-two nodes suppressed), or None if the family left no extant gene. to_newick serialises either.

origination is when the family was founded — the exact time of its origination event, or the root lineage's start for a family declared by initial_families. A :class:GeneNode records when it ended, so this is the one time the tree cannot derive: it is where the root's branch begins.

to_newick

to_newick(which: str = 'extant', *, annotate: bool = True) -> str | None

Newick of the "extant" (default) or "complete" tree; None if it is empty. Leaves are g<id>; with annotate internal nodes carry <kind>_n<species>; branch lengths are time differences.

The root carries one too, running from origination to where the root gene ended — the stem of the family, real time in which that founding gene existed. On the extant tree the root may be a node whose ancestors were suppressed; its branch still starts at origination and so absorbs them, exactly as the species tree's extant root absorbs its own.

Source code in zombi2/genomes/gene_trees.py

def to_newick(self, which: str = "extant", *, annotate: bool = True) -> str | None:
    """Newick of the ``"extant"`` (default) or ``"complete"`` tree; ``None`` if it is empty.
    Leaves are ``g<id>``; with ``annotate`` internal nodes carry ``<kind>_n<species>``; branch
    lengths are time differences.

    The root carries one too, running from ``origination`` to where the root gene ended — the
    stem of the family, real time in which that founding gene existed. On the extant tree the
    root may be a node whose ancestors were suppressed; its branch still starts at ``origination``
    and so absorbs them, exactly as the species tree's extant root absorbs its own."""
    root = self.extant if which == "extant" else self.complete
    if root is None:
        return None
    return _to_newick(root, annotate, self.origination) + ";"

zombi2.genomes.GeneCopy dataclass

GeneCopy(id: int, family: int)

One gene copy: a member of family family, identified by a globally-unique id. Its birth/death times and parentage live in the event log (the source of truth); the copy carries only what a genome snapshot needs to be self-describing — who it is and which family it is in. A genome may hold several copies sharing a family (that family's copy count).

Sequences

zombi2.sequences.simulate_sequences

simulate_sequences(genomes, *, model: SubstitutionModel, length: int | None = None, intergene_model: SubstitutionModel | None = None, intergene_speed=3.0, substitution=1.0, seed=None, parallel=False, progress=False) -> SequencesResult

Evolve one sequence down each family's gene tree under a substitution model.

genomes is a genome run — the :class:~zombi2.genomes.FamilyGenomesResult that genomes.simulate_genomes_family(...) returned. Its gene_trees are what the sequences evolve along and its complete_tree is the species tree the lineage clock rides; bare gene trees are rejected (they would run, but with no clock and no species phylogram — a silent degradation). Each family's complete gene tree is evolved, so the true history is complete and ancestral sequences exist for extinct/lost lineages too; the observable alignments are the extant tips.

model is a substitution model from the menu (:mod:.substitution_models) — nucleotide jc69 · k80 · hky85 · gtr, or protein poisson · jtt · dayhoff · wag · lg; its alphabet is what the sequences are written in (ACGT or the 20 amino acids). length is the number of sites. substitution is the per-site substitution rate (default 1.0): a branch of Δt time accrues substitution · Δt substitutions/site. The founding sequence of each family is drawn from the model's stationary frequencies. Deterministic given seed.

substitution may carry a lineage clock — one factor per species branch, shared across families, computed once before evolving, rescaling each gene-tree branch by the clock of the species branch it sits on: 1.0 * mod.ByLineage(spread=) is the uncorrelated clock (each branch drawn i.i.d.), and 1.0 * mod.FromParent(spread=) is the autocorrelated clock (the factor drifts parent→child down the species tree). Any other modifier (the Markov clock, the ByFamily per-family speed, +Γ) or a non-PerSite scope raises.

On a nucleotide genome run every root block is evolved — spacer as well as genes — each at its own length in bp, so length does not apply and is rejected. model evolves the genes and intergene_model (default jc69) the spacer, at intergene_speed times the rate (default 3.0). Because the whole genome is covered, the run also puts the genomes back together: .genomes holds every node's chromosomes, blocks concatenated in physical order — the complete tree, reconstructed — and .initial_genome the one the run started with.

The result carries the phylograms the sequences were drawn along — each gene tree and the species tree, with branch lengths converted from time to substitutions/site by the same base × clock × Δt.

parallel opts into evolving the gene trees concurrently — one gene tree per worker process, which is where a run's time goes when the sequences are long or a nucleotide genome has thousands of blocks. False (the default) runs the serial engine above. True uses every core; a positive int sets the worker count. It is a separate engine: each family draws from its own RNG stream (spawned from seed), so every worker count returns the same bytes, but that realisation differs from the serial one for a given seed — parallel is a speed choice, made once, not a drop-in for the default. Threads would not help here (numpy releases the GIL too little for these array sizes), so this is process-backed and only pays off above a work threshold. Because it spawns processes, a script that calls it with parallel set must guard its entry with if __name__ == "__main__": (the standard multiprocessing requirement); the zombi2 CLI already does, so --parallel there needs nothing extra.

Source code in zombi2/sequences/__init__.py

def simulate_sequences(genomes, *, model: SubstitutionModel, length: int | None = None,
                       intergene_model: SubstitutionModel | None = None, intergene_speed=3.0,
                       substitution=1.0, seed=None, parallel=False, progress=False) -> SequencesResult:
    """Evolve one sequence down each family's gene tree under a substitution ``model``.

    ``genomes`` is a **genome run** — the :class:`~zombi2.genomes.FamilyGenomesResult` that
    ``genomes.simulate_genomes_family(...)`` returned. Its ``gene_trees`` are what the sequences
    evolve along and its ``complete_tree`` is the species tree the lineage clock rides; bare gene
    trees are rejected (they would run, but with no clock and no species phylogram — a silent
    degradation). Each family's *complete* gene tree is evolved, so the true history is complete and
    ancestral sequences exist for extinct/lost lineages too; the observable ``alignments`` are the
    extant tips.

    ``model`` is a substitution model from the menu (:mod:`.substitution_models`) — nucleotide
    ``jc69`` · ``k80`` · ``hky85`` · ``gtr``, or protein ``poisson`` · ``jtt`` · ``dayhoff`` ·
    ``wag`` · ``lg``; its alphabet is what the sequences are written in (``ACGT`` or the 20 amino
    acids). ``length`` is the number of sites. ``substitution`` is the per-site substitution
    rate (default ``1.0``): a branch of ``Δt`` time accrues ``substitution · Δt`` substitutions/site.
    The founding sequence of each family is drawn from the model's stationary frequencies. Deterministic
    given ``seed``.

    ``substitution`` may carry a **lineage clock** — one factor per species branch, shared across
    families, computed once before evolving, rescaling each gene-tree branch by the clock of the species
    branch it sits on: ``1.0 * mod.ByLineage(spread=)`` is the uncorrelated clock (each branch drawn
    i.i.d.), and ``1.0 * mod.FromParent(spread=)`` is the autocorrelated clock (the factor drifts
    parent→child down the species tree). Any other modifier (the ``Markov`` clock, the ``ByFamily``
    per-family speed, ``+Γ``) or a non-``PerSite`` scope raises.

    On a **nucleotide** genome run every root block is evolved — spacer as well as genes — each at its
    own length in bp, so ``length`` does not apply and is rejected. ``model`` evolves the genes and
    ``intergene_model`` (default ``jc69``) the spacer, at ``intergene_speed`` times the rate (default
    ``3.0``). Because the whole genome is covered, the run also **puts the genomes back together**:
    ``.genomes`` holds every node's chromosomes, blocks concatenated in physical order — the complete
    tree, reconstructed — and ``.initial_genome`` the one the run started with.

    The result carries the **phylograms** the sequences were drawn along — each gene tree and the
    species tree, with branch lengths converted from time to substitutions/site by the same
    ``base × clock × Δt``.

    ``parallel`` opts into evolving the gene trees **concurrently** — one gene tree per worker
    process, which is where a run's time goes when the sequences are long or a nucleotide genome has
    thousands of blocks. ``False`` (the default) runs the serial engine above. ``True`` uses every
    core; a positive ``int`` sets the worker count. It is a **separate engine**: each family draws
    from its own RNG stream (spawned from ``seed``), so every worker count returns the *same* bytes,
    but that realisation differs from the serial one for a given seed — parallel is a speed choice,
    made once, not a drop-in for the default. Threads would not help here (numpy releases the GIL too
    little for these array sizes), so this is process-backed and only pays off above a work threshold.
    Because it spawns processes, a script that calls it with ``parallel`` set must guard its entry with
    ``if __name__ == "__main__":`` (the standard multiprocessing requirement); the ``zombi2`` CLI
    already does, so ``--parallel`` there needs nothing extra.
    """
    from ..genomes import NucleotideGenomesResult

    nucleotide = isinstance(genomes, NucleotideGenomesResult)
    if not nucleotide and not isinstance(genomes, FamilyGenomesResult):
        raise TypeError(
            f"the sequence level runs on a genome run, got {type(genomes).__name__} — pass the "
            "FamilyGenomesResult that genomes.simulate_genomes_family(...) returned, or the "
            "NucleotideGenomesResult from simulate_genomes_nucleotide(...): the whole run, not its "
            ".gene_trees. A sequence lives inside a gene, but its clock rides the *species* branch "
            "that gene sits on — one draw per lineage, shared by every family — so the run needs "
            "the species tree too."
        )
    species_tree = genomes.complete_tree
    if not isinstance(model, SubstitutionModel):
        raise TypeError(f"model must be a SubstitutionModel (e.g. hky85(kappa=2.0)), got {model!r}")
    if intergene_model is not None and not isinstance(intergene_model, SubstitutionModel):
        raise TypeError(f"intergene_model must be a SubstitutionModel, got {intergene_model!r}")

    if nucleotide:
        # Every recovered root block evolves — spacer as well as genes — so the run reconstructs the
        # whole genome rather than the declared loci. Each block brings its own length in bp, which
        # is why a single `length` would contradict the coordinates the genome recorded.
        if length is not None:
            raise ValueError(
                "length does not apply to a nucleotide genome run: every block carries its own "
                "length in bp, so one number here would contradict the coordinates the genomes run "
                "wrote. Drop it — the genome sets the lengths.")
        for name, m in (("model", model), ("intergene_model", intergene_model)):
            if m is not None and m.alphabet != BASES:
                raise ValueError(
                    f"{name}={m.name} is a protein model, but a nucleotide genome is measured in base "
                    "pairs and its blocks are read on either strand — amino acids have no complement "
                    "to read back. Use a nucleotide model (jc69 / k80 / hky85 / gtr).")
        if intergene_model is None:
            intergene_model = jc69()          # flat and parameterless: the null for unconstrained DNA
        if isinstance(intergene_speed, bool) or not isinstance(intergene_speed, (int, float)) \
                or intergene_speed <= 0:
            raise ValueError(f"intergene_speed must be a positive number, got {intergene_speed!r}")
        gene_trees = genomes.block_trees
        blocks = genomes.root_blocks
        genic = {span: fam for fam, span in genomes.gene_spans.items()}
        # per block: its length, whether it is genic, the model it evolves under and its speed
        per_block = {}
        for i, (src, a, b) in enumerate(blocks):
            is_gene = (src, a, b) in genic
            per_block[i] = (b - a, model if is_gene else intergene_model,
                            1.0 if is_gene else float(intergene_speed))
        # Founded from a real FASTA: a block's founding sequence is the supplied DNA at its own root
        # coordinates, encoded to states, rather than a stationary draw. A de-novo source is not in
        # initial_sequence (it arose mid-run), so its blocks still draw from the model. `None` per block
        # ⇒ draw, exactly as before, so a run without one is unchanged.
        founding_seed: dict[int, "np.ndarray | None"] = {}
        for i, (src, a, b) in enumerate(blocks):
            root = genomes.initial_sequence.get(src)
            if root is None:
                founding_seed[i] = None
                continue
            f_model = model if (src, a, b) in genic else intergene_model
            if f_model.alphabet != BASES:
                raise ValueError(
                    f"the run was founded from a FASTA (DNA), but {f_model.name} is a protein model — "
                    "a nucleotide sequence cannot found an amino-acid alignment")
            founding_seed[i] = encode(root[a:b], f_model.alphabet)
    else:
        gene_trees = genomes.gene_trees
        if length is None:
            raise ValueError("length is required: the number of sites each family evolves")
        if isinstance(length, bool) or not isinstance(length, int) or length < 1:
            raise ValueError(f"length must be a positive integer, got {length!r}")
        if intergene_model is not None:
            raise ValueError(
                "intergene_model applies to a nucleotide genome run, where blocks are genes or "
                "spacer. A family or ordered run has gene families only, so there is nothing "
                "for a second model to evolve.")
        per_block = None
    rate = as_rate(substitution, default_scope=PerSite)
    if not isinstance(rate.scope, PerSite):
        raise ValueError(
            f"substitution has a {type(rate.scope).__name__} scope, but the sequence engine wires only "
            f"PerSite (the default) this slice — drop the scope wrapper or use PerSite(...)."
        )
    clock_mod = None
    if rate.modifiers:
        if len(rate.modifiers) == 1 and isinstance(rate.modifiers[0], (ByLineage, FromParent)):
            clock_mod = rate.modifiers[0]
        else:
            offenders = ", ".join(sorted({type(m).__name__ for m in rate.modifiers
                                          if not isinstance(m, (ByLineage, FromParent))})
                                  or ["a second clock"])
            raise ValueError(
                f"substitution carries {offenders}, but this slice wires a single lineage clock — one "
                "ByLineage (uncorrelated) or one FromParent (autocorrelated). The Markov clock, the "
                "ByFamily per-family speed, and +Γ across-site heterogeneity are not wired."
            )
    rate_base = rate.base

    alignments: dict[int, dict[str, str]] = {}
    ancestral: dict[int, dict[str, str]] = {}
    founding: dict[int, str] = {}
    phylograms: dict[int, dict[str, str | None]] = {}
    if not parallel:
        # Serial reference engine — the default, left exactly as it was. One shared generator draws the
        # clock, then each family is walked in turn. `parallel` selects a *separate* engine (decision A),
        # so turning it on gives a different-but-valid realisation for a seed; this path never changes.
        rng = np.random.default_rng(seed)
        clock = _draw_clock(clock_mod, species_tree, gene_trees, rng)
        # One transition-CDF cache per model, shared across every block that model evolves. Branch lengths
        # recur across blocks (a block passing straight through a species branch reuses its length), so a
        # run-wide cache builds a few hundred matrices where a per-block cache rebuilt tens of thousands.
        # Keyed by model identity — genes and spacer are different models and must not share a cache.
        cdf_caches: dict[int, dict[float, np.ndarray]] = {}
        bar = progress_bar(len(gene_trees), "sequences", unit="family", enabled=progress)
        for family in sorted(gene_trees):  # sorted for reproducibility given the seed
            bar.update()
            gt = gene_trees[family]
            if per_block is None:
                f_len, f_model, f_rate = length, model, rate_base
            else:                       # a nucleotide block: its own length, and spacer runs faster
                f_len, f_model, speed = per_block[family]
                f_rate = rate_base * speed
            seed_states = None if per_block is None else founding_seed[family]
            cache = cdf_caches.setdefault(id(f_model), {})
            states, founding_states = evolve_gene_tree(gt.complete, f_model, f_len, f_rate, clock, rng,
                                                       gt.origination, founding=seed_states,
                                                       cdf_cache=cache)
            alignments[family], ancestral[family] = _split(gt, states, f_model)
            founding[family] = decode(founding_states, f_model.alphabet)
            scaled = _scaled_gene_tree(gt, f_rate, clock)  # branch lengths in subs/site
            ext = scaled.extant
            phylograms[family] = {"complete": _gene_newick(scaled.complete),
                                  "extant": _gene_newick(ext) if ext is not None else None}
        bar.close()
    else:
        # Parallel engine (opt-in): one gene tree per process, each under its own spawned RNG stream, so
        # any worker count is bit-identical to any other. The clock is a shared read-only draw, so it is
        # taken once here from a reserved stream (index 0) and shipped to every worker.
        from .._runtime.parallel import guard_pool_workers, resolve_workers
        from ._pergenetree import evolve_families
        workers = guard_pool_workers(resolve_workers(parallel))
        spawned = np.random.SeedSequence(seed).spawn(1 + len(gene_trees))
        clock = _draw_clock(clock_mod, species_tree, gene_trees, np.random.default_rng(spawned[0]))
        alignments, ancestral, founding, phylograms = evolve_families(
            gene_trees, per_block, model, intergene_model, length, rate_base, clock,
            founding_seed if nucleotide else None, spawned[1:], workers, progress)

    sp_scaled = _scaled_species_tree(species_tree, rate_base, clock)   # the clock made visible
    sp_extant = prune(sp_scaled, keep="extant")
    species_phylogram = {"complete": sp_scaled.to_newick(),
                         "extant": sp_extant.to_newick() if sp_extant is not None else None}
    # A nucleotide run evolved every block, so **every** node's genome can be put back together —
    # one map, as at the genome level. Which sequences each node reads is the split the level already
    # makes: an extant tip's genes are tips of their block trees, everything else's are not. The
    # concatenation is deferred to read-time (see :class:`_AssembledGenomes`): only the cheap per-node
    # layout is captured now, so hundreds of megabases of genome do not all sit in memory beside the
    # per-block sequences they are built from.
    assembled: "Mapping[str, dict[int, str]]" = {}
    initial_genome: dict[int, str] = {}
    if nucleotide:
        # Capture the layouts in the same order the eager build used — extant nodes (read from
        # `alignments`) sorted first, then the rest (read from `ancestral`) — so the map iterates,
        # and `write` emits its files, in exactly the previous order.
        extant_ids = sorted(n.id for n in species_tree.extant())
        extant_id_set = set(extant_ids)
        extant_labels = {node_label(i) for i in extant_ids}
        ordered_ids = extant_ids + [i for i in sorted(species_tree.nodes) if i not in extant_id_set]
        layouts = {node_label(i): genomes.assembly(i) for i in ordered_ids}
        assembled = _AssembledGenomes(layouts, alignments, ancestral, extant_labels)
        # The genome the run started with. Its blocks were all laid down at the start, so each one's
        # sequence there is its `founding` draw — the state the stem leads *from*. It is not a node,
        # so it is in neither map above; the same reason `founding` is not in `ancestral`.
        for cid, pieces in genomes.initial_assembly().items():
            initial_genome[cid] = "".join(
                founding[block] if strand == 1 else founding[block].translate(_COMPLEMENT)[::-1]
                for (block, strand) in pieces)

    return SequencesResult(alignments, ancestral, founding, phylograms, species_phylogram, seed,
                           assembled, initial_genome, "block" if nucleotide else "family")

zombi2.sequences.SequencesResult dataclass

SequencesResult(alignments: dict[int, dict[str, str]], ancestral: dict[int, dict[str, str]], founding: dict[int, str], phylograms: dict[int, dict[str, str | None]], species_phylogram: dict[str, str | None], seed: int | None, genomes: 'Mapping[str, dict[int, str]]' = dict(), initial_genome: dict[int, str] = dict(), unit: str = 'family')

What :func:simulate_sequences returns.

• alignments — {family: {g<copy>: sequence}}: the observable gene alignment, one entry per extant gene-tree tip, keyed by its (unique, per-segment) gene id — the same labels as the gene tree's / phylogram's Newick leaves. Empty for a family with no surviving copy.
• ancestral — {family: {g<copy>: sequence}}: the true sequence at every node that is not an extant tip — internal nodes (the family's root gene included) and the dead tips, where a copy was lost or its species went extinct. With alignments it accounts for every node of the tree exactly once, so every label in the complete phylogram names a sequence.
• founding — {family: sequence}: the sequence the family started with, at its origination — the state the phylogram's root branch leads from. It is drawn from the model's stationary frequencies and then evolves across the stem into the root gene's sequence, so it is not the same string as ancestral[family]["g<root copy>"] unless the stem is empty. Kept out of ancestral on purpose: those keys pair one-to-one with phylogram nodes, and the origination is a point on a branch, not a node.
• phylograms — {family: {"complete": newick, "extant": newick | None}}: each gene tree with branch lengths in substitutions/site (base × lineage-clock × Δt) — the ground-truth tree behind each alignment. Every node is labelled by its gene id g<copy>, so the tips match the alignments keys and the internal nodes match the ancestral keys (the phylogram pairs one-to-one with the sequences). "extant" is None for a family with no survivor.
• species_phylogram — {"complete": newick, "extant": newick | None}: the species tree with branch lengths in substitutions/site — the molecular clock made visible (which lineages ran hot / cold). Always present: a run always comes from a genome run, which carries its tree.
• genomes — {lineage: {chromosome id: sequence}}: every node's assembled genome, its blocks concatenated in physical order (reverse-complemented where the genome carries them inverted) — extant tips, ancestors and the lineages that went extinct alike. The same coverage as the genome level's own genomes, which is keyed the same way: the observed ones are the extant tips, {node_label(n.id) for n in complete_tree.extant()}. Only a nucleotide genome run has any — a family or ordered run has gene families, not coordinates, so there is no genome to lay out and this is empty.
• initial_genome — {chromosome id: sequence}: the genome the run started with, at the root lineage's origination. Not in genomes, because it belongs to no node: the root branch is real simulated time, so the root node's genome is this one plus whatever happened along the stem. It stands to genomes as founding stands to ancestral.
• seed — the run's seed.
• unit — what the integer key of alignments / ancestral / founding / phylograms names: "family" (a gene family id) on a family or ordered run, "block" (an index into the genome run's root_blocks) on a nucleotide one, where every block evolves and spacer has no family. They are different numbering schemes over the same ints, so a gene family id is not a key here on a nucleotide run — go through :meth:~zombi2.genomes.NucleotideGenomesResult.block_of. It is also what the filenames say.

write

write(directory, outputs=('alignments', 'phylograms', 'species_phylogram', 'genomes', 'initial_genome')) -> None

Write chosen outputs to directory (created if needed). <u> below is fam<family> on a family or ordered run and block<index> on a nucleotide one — the integer keys mean different things, so the files say which (see :attr:unit):

• "alignments" → <u>.fasta (skipped for empty families).
• "ancestral" → sequences_ancestral_<u>.fasta.
• "founding" → sequences_founding.fasta, one record <u> apiece: the sequence each family originated with, before its stem.
• "phylograms" → phylogram_<u>_{complete,extant}.nwk (subs/site).
• "species_phylogram" → clock_species_tree_{complete,extant}.nwk: the species tree with its branches in substitutions/site — the molecular clock made visible.
• "genomes" → genome_<lineage>.fasta, one file per node — extant, extinct and ancestral alike — with one record per chromosome. Nucleotide runs only; nothing is written otherwise. The big one: a real genome times every node in the tree.
• "initial_genome" → genome_initial.fasta, the genome the run started with.

Source code in zombi2/sequences/__init__.py

def write(self, directory,
          outputs=("alignments", "phylograms", "species_phylogram", "genomes",
                   "initial_genome")) -> None:
    """Write chosen ``outputs`` to ``directory`` (created if needed). ``<u>`` below is
    ``fam<family>`` on a family or ordered run and ``block<index>`` on a nucleotide one — the
    integer keys mean different things, so the files say which (see :attr:`unit`):

    - ``"alignments"`` → ``<u>.fasta`` (skipped for empty families).
    - ``"ancestral"`` → ``sequences_ancestral_<u>.fasta``.
    - ``"founding"`` → ``sequences_founding.fasta``, one record ``<u>`` apiece: the sequence each
      family originated with, before its stem.
    - ``"phylograms"`` → ``phylogram_<u>_{complete,extant}.nwk`` (subs/site).
    - ``"species_phylogram"`` → ``clock_species_tree_{complete,extant}.nwk``: the species tree
      with its branches in substitutions/site — the molecular clock made visible.
    - ``"genomes"`` → ``genome_<lineage>.fasta``, one file per node — extant, extinct and
      ancestral alike — with one record per chromosome. Nucleotide runs only; nothing is written
      otherwise. The big one: a real genome times every node in the tree.
    - ``"initial_genome"`` → ``genome_initial.fasta``, the genome the run started with.
    """
    unknown = [o for o in outputs if o not in _WRITE_OUTPUTS]
    if unknown:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(_WRITE_OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    u = self._stem
    if "alignments" in outputs:
        for fam, aln in self.alignments.items():
            if aln:
                _write_fasta(d / f"{u}{fam}.fasta", aln)
    if "ancestral" in outputs:
        for fam, anc in self.ancestral.items():
            if anc:
                _write_fasta(d / f"sequences_ancestral_{u}{fam}.fasta", anc)
    if "founding" in outputs and self.founding:
        _write_fasta(d / "sequences_founding.fasta",
                     {f"{u}{fam}": seq for fam, seq in sorted(self.founding.items())})
    if "phylograms" in outputs:
        for fam, ph in self.phylograms.items():
            (d / f"phylogram_{u}{fam}_complete.nwk").write_text(ph["complete"] + "\n")
            if ph["extant"] is not None:
                (d / f"phylogram_{u}{fam}_extant.nwk").write_text(ph["extant"] + "\n")
    if "species_phylogram" in outputs:
        sp = self.species_phylogram
        (d / "clock_species_tree_complete.nwk").write_text(sp["complete"] + "\n")
        if sp["extant"] is not None:
            (d / "clock_species_tree_extant.nwk").write_text(sp["extant"] + "\n")
    # every genome is written the same way and named by whose it is — a node label, or "initial"
    for token, genomes in (("genomes", self.genomes),
                           ("initial_genome",
                            {"initial": self.initial_genome} if self.initial_genome else {})):
        if token in outputs:
            for lineage, chroms in genomes.items():
                _write_fasta(d / f"genome_{lineage}.fasta",
                             {f"{lineage}_chr{cid}": seq for cid, seq in chroms.items()})

Traits

zombi2.traits.simulate_continuous

simulate_continuous(tree, *, start=0.0, rate=1.0, reverts_to=None, pull=None, correlation=None, at_speciation=None, regimes=None, seed=None, progress=False) -> TraitsResult

Evolve a continuous trait down a tree and return a :class:TraitsResult. One process, its variants selected by knobs (SPEC §4): Brownian motion (bare rate), Ornstein–Uhlenbeck (add reverts_to + pull), early burst (a OnTime skyline on rate), and variable-rates BM (an FromParent modifier on rate).

Correlated traits ride together in one call (the joint rule inside a level): pass start and rate as dicts keyed by trait name and a correlation={(a, b): ρ} overlay (each ρ ∈ [−1, 1]). The traits then diffuse jointly — the branch increment is drawn from MVN(0, Σ·dt) with Σ = D R D (D = diag(σ_i), R the correlation matrix), so at a tip the correlation between two traits is exactly their ρ. Correlated Brownian motion only, with bare per-trait rates.

tree is the complete species tree (a :class:~zombi2.species.Tree, or a :class:~zombi2.species.SpeciesResult whose complete_tree is used). The trait evolves on every lineage, extant and extinct alike, so the ancestral states are exact and complete; the observed dataset is the extant tips, result.values.

start is the value at t = 0 (the crown origin, root.birth_time): the root lineage diffuses over its own branch [0, first split] like any other, so a trait and a genome evolve over the same branch set, and each node's stored value is the trait at that node's end_time (node_values[root] is the value at the first split, not start).

rate is the variance-rate σ² (a scope(base) × modifiers rate spec), per lineage: each lineage diffuses independently at σ², never pooled across the tree. A bare number is Brownian motion (Normal(0, σ²·dt) over a branch); a OnTime modifier makes σ² change through time — early burst / ACDC — with the per-branch variance the exact integral ∫ σ²(t) dt; an FromParent(spread=…) modifier makes σ² drift branch-to-branch — variable-rates BM ("ClaDS for traits") — each lineage inheriting its parent's σ² times a lognormal kick drawn at the split; a OnTotalDiversity(cap=…) modifier makes σ² slow as the clade fills up — diversity-dependent / ecological-limits trait evolution — σ² scaled by (1 − standing_diversity/cap) as the tree's lineages-through-time grows (the tree is a fixed input the trait reads).

reverts_to (the optimum θ) and pull (the strength α > 0) turn the diffusion into Ornstein–Uhlenbeck — the value is pulled toward θ while it diffuses, the exact per-branch transition being Normal(θ + (x−θ)·e^{−α·dt}, σ²/(2α)·(1−e^{−2α·dt})). Give both or neither. OU with a modified σ² (early burst or variable rates on rate) is not wired yet — use one or the other.

at_speciation adds an on-speciation jump — Normal(0, at_speciation) on each daughter at every speciation (the punctuational mode), layered on top of the along-branch anagenesis. regimes gives multi-optimum OU: pass a discrete :class:TraitsResult (a stochastic map painted by :func:simulate_discrete on this same tree) and a per-regime reverts_to={regime: θ}, and the value follows OU toward whichever regime's optimum a branch is in. Deterministic given seed.

Source code in zombi2/traits/continuous.py

def simulate_continuous(tree, *, start=0.0, rate=1.0, reverts_to=None, pull=None,
                        correlation=None, at_speciation=None, regimes=None, seed=None,
                        progress=False) -> TraitsResult:
    """Evolve a continuous trait down a tree and return a :class:`TraitsResult`. One process, its
    variants selected by knobs (SPEC §4): **Brownian motion** (bare ``rate``), **Ornstein–Uhlenbeck**
    (add ``reverts_to`` + ``pull``), **early burst** (a ``OnTime`` skyline on ``rate``), and
    **variable-rates BM** (an ``FromParent`` modifier on ``rate``).

    **Correlated traits** ride together in **one call** (the joint rule inside a level): pass
    ``start`` and ``rate`` as dicts keyed by trait name and a ``correlation={(a, b): ρ}`` overlay
    (each ρ ∈ [−1, 1]). The traits then diffuse jointly — the branch increment is drawn from
    ``MVN(0, Σ·dt)`` with ``Σ = D R D`` (``D = diag(σ_i)``, ``R`` the correlation matrix), so at a tip
    the correlation between two traits is exactly their ρ. Correlated **Brownian motion** only, with bare per-trait rates.

    ``tree`` is the **complete** species tree (a :class:`~zombi2.species.Tree`, or a
    :class:`~zombi2.species.SpeciesResult` whose ``complete_tree`` is used). The trait evolves on
    **every** lineage, extant and extinct alike, so the ancestral states are exact and complete; the
    observed dataset is the extant tips, ``result.values``.

    ``start`` is the value at ``t = 0`` (the crown origin, ``root.birth_time``): the root lineage
    diffuses over its own branch ``[0, first split]`` like any other, so a trait and a genome evolve
    over the **same** branch set, and each node's stored value is the trait at that node's
    ``end_time`` (``node_values[root]`` is the value at the first split, not ``start``).

    ``rate`` is the variance-rate σ² (a ``scope(base) × modifiers`` rate spec), *per lineage*: each
    lineage diffuses independently at σ², never pooled across the tree. A bare number is Brownian
    motion (``Normal(0, σ²·dt)`` over a branch); a ``OnTime`` modifier makes σ² change through time —
    early burst / ACDC — with the per-branch variance the exact integral ``∫ σ²(t) dt``; an
    ``FromParent(spread=…)`` modifier makes σ² **drift branch-to-branch** — variable-rates BM ("ClaDS
    for traits") — each lineage inheriting its parent's σ² times a lognormal kick drawn at the split;
    a ``OnTotalDiversity(cap=…)`` modifier makes σ² **slow as the clade fills up** — diversity-dependent /
    ecological-limits trait evolution — σ² scaled by ``(1 − standing_diversity/cap)`` as the tree's
    lineages-through-time grows (the tree is a fixed input the trait reads).

    ``reverts_to`` (the optimum θ) and ``pull`` (the strength α > 0) turn the diffusion into
    Ornstein–Uhlenbeck — the value is pulled toward θ while it diffuses, the exact per-branch
    transition being ``Normal(θ + (x−θ)·e^{−α·dt}, σ²/(2α)·(1−e^{−2α·dt}))``. Give **both** or
    neither. OU with a *modified* σ² (early burst or variable rates on ``rate``) is not wired yet —
    use one or the other.

    ``at_speciation`` adds an **on-speciation** jump — ``Normal(0, at_speciation)`` on each daughter at
    every speciation (the punctuational mode), layered on top of the along-branch anagenesis.
    ``regimes`` gives **multi-optimum OU**: pass a discrete :class:`TraitsResult` (a stochastic map
    painted by :func:`simulate_discrete` on this same tree) and a per-regime ``reverts_to={regime: θ}``,
    and the value follows OU toward whichever regime's optimum a branch is in. Deterministic given
    ``seed``.
    """
    tree = tree.complete_tree if isinstance(tree, SpeciesResult) else tree
    if regimes is not None:
        return _simulate_regimes(tree, start, rate, reverts_to, pull, regimes, at_speciation, seed,
                                 progress)
    if isinstance(start, dict) or isinstance(rate, dict) or correlation is not None:
        return _simulate_correlated(tree, start, rate, reverts_to, pull, correlation,
                                    at_speciation, seed, progress)
    if isinstance(start, bool) or not isinstance(start, (int, float)) or not math.isfinite(start):
        raise ValueError(f"start must be a finite number, got {start!r}")
    r = as_rate(rate, default_scope=PerLineage)
    if not isinstance(r.scope, PerLineage):
        raise ValueError(
            f"rate has a {type(r.scope).__name__} scope, but a continuous trait's variance-rate is "
            f"per lineage — drop the scope wrapper (per lineage is the default)."
        )
    # OnTime (early burst), FromParent (variable-rates BM), and OnTotalDiversity (diversity-dependent) are the
    # wired σ² modifiers; anything else is rejected loudly — the genome engine's discipline.
    for m in r.modifiers:
        if not isinstance(m, WIRED_MODIFIERS):
            raise ValueError(
                f"rate carries {type(m).__name__}, which the continuous trait engine does not support "
                f"— OnTime (early burst), FromParent (variable-rates BM), and OnTotalDiversity "
                f"(diversity-dependent) are wired."
            )
    drifts = [m for m in r.modifiers if isinstance(m, FromParent)]
    if len(drifts) > 1:
        raise ValueError("rate carries more than one FromParent modifier; a variance-rate drifts one way")
    drift = drifts[0] if drifts else None  # the per-lineage σ² drift (variable-rates BM), or None
    has_diversity = any(isinstance(m, OnTotalDiversity) for m in r.modifiers)  # σ² reads the standing LTT

    # OU: reverts_to (θ) + pull (α) turn the diffusion into mean-reversion — both or neither.
    is_ou = reverts_to is not None or pull is not None
    if is_ou:
        if reverts_to is None or pull is None:
            raise ValueError(
                "Ornstein–Uhlenbeck needs both reverts_to (the optimum) and pull (the strength); "
                "give both, or neither for Brownian motion."
            )
        if isinstance(reverts_to, bool) or not isinstance(reverts_to, (int, float)) \
                or not math.isfinite(reverts_to):
            raise ValueError(f"reverts_to must be a finite number, got {reverts_to!r}")
        if isinstance(pull, bool) or not isinstance(pull, (int, float)) \
                or not math.isfinite(pull) or pull <= 0:
            raise ValueError(
                f"pull must be a finite positive number (omit it for Brownian motion), got {pull!r}"
            )
        if r.modifiers:
            raise ValueError(
                "a modified variance-rate (early burst via OnTime, variable rates via FromParent, or "
                "diversity-dependence via OnTotalDiversity) combined with OU (reverts_to / pull) is not "
                "wired yet — use one or the other."
            )
        theta, alpha = float(reverts_to), float(pull)
        sigma2 = r.effective(lineages=1)  # σ² is constant under OU (modifiers are rejected above)

    jump_sd = _at_speciation_jump_sd(at_speciation)  # on-speciation jump width (0 if not requested)

    rng = np.random.default_rng(seed)
    ltt = _LTT(tree) if has_diversity else None  # the standing-diversity curve, when σ² reads it
    node_values: dict[int, float] = {}
    root = tree.nodes[tree.root]
    # the initial value at t=0 — the origin the log reconstructs from (SPEC §2). A diffusion cannot be
    # rebuilt from events, but the row keeps the file's shape uniform across trait kinds.
    events: list[Change] = [Change(root.birth_time, "root", tree.root, None, float(start))]
    inh: dict[int, float] = {}  # each lineage's σ² drift factor (variable-rates BM), constant per branch
    for i in _preorder(tree, progress):
        node = tree.nodes[i]
        # the root starts from `start` at t=0; every other node from its parent's end value (parent
        # < i, already set). One uniform rule: node_values[i] is the trait at node i's end_time.
        x = float(start) if node.parent is None else node_values[node.parent]
        if node.parent is not None and jump_sd > 0.0:
            jumped = x + float(rng.normal(0.0, jump_sd))  # on speciation: a jump at the split…
            events.append(Change(node.birth_time, "on_speciation", i, x, jumped))
            x = jumped                                    # …then anagenesis along the branch
        # thread the inherited factor: the root's is 1.0, each daughter's is its parent's times a
        # lognormal kick drawn at the split (so σ² is autocorrelated down the tree). None ⇒ 1.0, no draw.
        if node.parent is None:
            inh[i] = drift.initial() if drift else 1.0
        else:
            inh[i] = drift.descend(inh[node.parent], rng) if drift else 1.0
        t0, t1 = node.birth_time, node.end_time
        if is_ou:
            e = math.exp(-alpha * (t1 - t0))       # mean-reversion toward θ over the branch
            mean = theta + (x - theta) * e
            var = sigma2 / (2.0 * alpha) * (1.0 - e * e)
        else:
            mean = x                                # pure diffusion (BM / early burst / variable-rates)
            var = _accrued_variance(r, t0, t1, inherited=inh[i], ltt=ltt)
        std = math.sqrt(var) if var > 0.0 else 0.0
        node_values[i] = mean + (float(rng.normal(0.0, std)) if std > 0.0 else 0.0)

    return TraitsResult(tree, node_values, events, seed)

zombi2.traits.simulate_discrete

simulate_discrete(tree, *, states, switch=None, start=None, liability=None, threshold=None, correlation=None, at_speciation=None, seed=None, progress=False) -> TraitsResult

Evolve a discrete-state trait down a tree and return a :class:TraitsResult. Two mechanisms:

• Mk (switch=) — a continuous-time Markov chain over the states, simulated exactly by Gillespie along every branch, so each node's (state, duration) segments are the realized history (.history) and .events reads off the transitions. switch is a symmetric rate (0.1), a {"marine->terrestrial": 0.1} dict, or a k×k matrix (see :func:_q_matrix). start is the root state (a label in states; None draws one uniformly).
• Threshold (liability= + threshold=) — the Wright–Felsenstein model: a discrete state read off an underlying continuous Brownian liability (variance-rate liability), cut into states by the threshold cut point(s) (k−1 increasing cuts for k states). start is the initial liability (a number, default 0.0). Give liability as a dict + a correlation={(a, b): ρ} overlay to evolve correlated discrete traits jointly — their liabilities diffuse together (Σ = D R D) and each is cut by the shared thresholds. A threshold trait has no Gillespie map, so .history is None and .events empty.

tree is the complete species tree (a :class:~zombi2.species.Tree or :class:~zombi2.species.SpeciesResult); the trait evolves on every lineage (convention B: the root diffuses over its own branch), and .values reads the extant tips. On an Mk trait, at_speciation (a probability in [0, 1]) adds an on-speciation shift — each daughter hops to a uniformly-chosen other state with that chance at every speciation. Deterministic given seed.

Source code in zombi2/traits/discrete.py

def simulate_discrete(tree, *, states, switch=None, start=None, liability=None, threshold=None,
                      correlation=None, at_speciation=None, seed=None,
                      progress=False) -> TraitsResult:
    """Evolve a discrete-state trait down a tree and return a :class:`TraitsResult`. Two mechanisms:

    - **Mk** (``switch=``) — a continuous-time Markov chain over the ``states``, simulated **exactly**
      by Gillespie along every branch, so each node's ``(state, duration)`` segments *are* the realized
      history (``.history``) and ``.events`` reads off the transitions. ``switch`` is a symmetric rate
      (``0.1``), a ``{"marine->terrestrial": 0.1}`` dict, or a ``k×k`` matrix (see :func:`_q_matrix`).
      ``start`` is the root state (a label in ``states``; ``None`` draws one uniformly).
    - **Threshold** (``liability=`` + ``threshold=``) — the Wright–Felsenstein model: a discrete state
      read off an underlying continuous Brownian **liability** (variance-rate ``liability``), cut into
      ``states`` by the ``threshold`` cut point(s) (``k−1`` increasing cuts for ``k`` states). ``start``
      is the initial *liability* (a number, default 0.0). Give ``liability`` as a dict + a
      ``correlation={(a, b): ρ}`` overlay to evolve **correlated** discrete traits jointly — their
      liabilities diffuse together (``Σ = D R D``) and each is cut by the shared thresholds. A threshold
      trait has no Gillespie map, so ``.history`` is ``None`` and ``.events`` empty.

    ``tree`` is the **complete** species tree (a :class:`~zombi2.species.Tree` or
    :class:`~zombi2.species.SpeciesResult`); the trait evolves on every lineage (convention B: the root
    diffuses over its own branch), and ``.values`` reads the extant tips. On an Mk trait,
    ``at_speciation`` (a probability in ``[0, 1]``) adds an **on-speciation** shift — each daughter hops
    to a uniformly-chosen other state with that chance at every speciation. Deterministic given ``seed``.
    """
    tree = tree.complete_tree if isinstance(tree, SpeciesResult) else tree
    states = list(states)
    if len(states) < 2:
        raise ValueError(f"a discrete trait needs at least 2 states, got {states!r}")
    if len(set(states)) != len(states):
        raise ValueError(f"states must be unique, got {states!r}")
    if liability is not None or threshold is not None:
        if switch is not None:
            raise ValueError("give switch= (an Mk trait) OR liability=/threshold= (a threshold trait), not both")
        if at_speciation is not None:
            raise ValueError("at_speciation is not wired for threshold traits yet — it applies to Mk (switch=) traits")
        return _simulate_threshold(tree, states, liability, threshold, start, correlation, seed,
                                   progress)
    if correlation is not None:
        raise ValueError("correlation= on a discrete trait needs the threshold model — give liability= and threshold=")
    if switch is None:
        raise ValueError("give switch= — the transition rate(s) between the discrete states.")
    Q = _q_matrix(states, switch)
    if at_speciation is not None and (isinstance(at_speciation, bool)
            or not isinstance(at_speciation, (int, float)) or not 0.0 <= at_speciation <= 1.0):
        raise ValueError(f"at_speciation must be a probability in [0, 1] (the shift chance), got {at_speciation!r}")
    shift = 0.0 if at_speciation is None else float(at_speciation)

    rng = np.random.default_rng(seed)
    idx = {s: i for i, s in enumerate(states)}
    if start is None:
        start_i = int(rng.integers(len(states)))
    elif start in idx:
        start_i = idx[start]
    else:
        raise ValueError(
            f"start must be one of states={states} (or None for a uniform draw), got {start!r}"
        )

    node_values: dict[int, object] = {}
    root = tree.nodes[tree.root]
    # the initial state at t=0 — the origin the log reconstructs from: tree + this + the switches give
    # the driver on every lineage, so the event log is the conditioning file (no separate driver).
    events: list[Change] = [Change(root.birth_time, "root", tree.root, None, states[start_i])]
    for i in _preorder(tree, progress):
        node = tree.nodes[i]
        # the root starts from `start` at t=0 and evolves over its own branch; every other node from
        # its parent's end state (parent < i, already set) — the same convention-B walk as continuous.
        cur = start_i if node.parent is None else idx[node_values[node.parent]]
        if node.parent is not None and shift > 0.0 and float(rng.random()) < shift:
            j = int(rng.integers(len(states) - 1))  # on speciation: hop to a uniform *other* state
            new = j if j < cur else j + 1
            events.append(Change(node.birth_time, "on_speciation", i, states[cur], states[new]))
            cur = new
        end_i, segs = _gillespie(cur, node.end_time - node.birth_time, Q, rng)
        t = node.birth_time  # the transitions between the Gillespie segments are the on-branch events
        for (s1, d1), (s2, _d) in zip(segs, segs[1:]):
            t += d1
            events.append(Change(t, "on_branch", i, states[s1], states[s2]))
        node_values[i] = states[end_i]
    events.sort(key=lambda c: c.time)
    return TraitsResult(tree, node_values, events, seed, kind="discrete")

zombi2.traits.TraitsResult dataclass

TraitsResult(complete_tree: Tree, node_values: dict[int, object], events: list[Change] = list(), seed: int | None = None, kind: str = 'continuous')

What simulate_continuous / simulate_discrete returns: the complete_tree it ran on, node_values at every node (the value at each node — extant, extinct, and internal alike; a float for a continuous trait, a state label for a discrete / threshold one, a per-trait dict for correlated traits), the events log, the seed, and the kind ("continuous" / "discrete" / "threshold"). The observed dataset is the extant tips, .values.

events is the timestamped event log — the same shape as the genome level's and the source of truth for a discrete (Mk) trait, from which history (the per-branch stochastic character map) is derived. A continuous trait has no along-branch events, so its log holds only the on-speciation jumps (empty without at_speciation) while node_values carries the diffusion; a threshold trait's crossings are un-timed, so its log is empty and it has no map.

values property

values: dict[int, object]

The observed trait dataset — the value at each extant tip (the comparative-data vector). Internal and extinct nodes keep their exact ancestral / lineage values in node_values.

history cached property

history: dict[int, list] | None

The per-branch stochastic character map — {node: [(state, duration), …]} whose durations sum to the branch length — derived from the event log (a discrete / Mk trait only). None for a continuous trait (a diffusion has no map) and for a threshold trait (its liability crossings are un-timed).

write

write(directory, outputs=('values',)) -> None

Write chosen outputs to directory (created if needed): "values" → trait_values.tsv (the node<TAB>trait table over the extant tips); "events" → trait_events.tsv, the event log (time · kind · lineage · from · to) — one root row at the crown giving the initial state, then every switch in time order; "tree" → trait_tree.nwk, the complete tree as Newick with every node annotated [&trait=…] (a trait tree, carrying the exact ancestral values; opens in FigTree / iTOL).

trait_events.tsv is also the conditioning file: a genome / sequence run drives a rate with mod.DrivenBy("trait_events.tsv", …), replaying it against the shared tree. A discrete trait's log reconstructs its state on every lineage exactly (that is what the root row and the switch times are for); a continuous trait's diffusion cannot be rebuilt from events, so it carries only the root row and any on-speciation jumps.

Source code in zombi2/traits/result.py

def write(self, directory, outputs=("values",)) -> None:
    """Write chosen ``outputs`` to ``directory`` (created if needed): ``"values"`` →
    ``trait_values.tsv`` (the ``node<TAB>trait`` table over the extant tips); ``"events"`` →
    ``trait_events.tsv``, the event log (``time · kind · lineage · from · to``) — one ``root`` row
    at the crown giving the initial state, then every switch in time order; ``"tree"`` →
    ``trait_tree.nwk``, the complete tree as Newick with **every** node annotated ``[&trait=…]``
    (a *trait tree*, carrying the exact ancestral values; opens in FigTree / iTOL).

    ``trait_events.tsv`` is also the **conditioning file**: a genome / sequence run drives a rate
    with ``mod.DrivenBy("trait_events.tsv", …)``, replaying it against the shared tree. A
    **discrete** trait's log reconstructs its state on every lineage exactly (that is what the
    ``root`` row and the switch times are for); a continuous trait's diffusion cannot be rebuilt
    from events, so it carries only the ``root`` row and any on-speciation jumps."""
    unknown = [o for o in outputs if o not in _WRITE_OUTPUTS]
    if unknown:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(_WRITE_OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    if "values" in outputs:
        (d / "trait_values.tsv").write_text(_values_tsv(self.values))
    if "events" in outputs:
        (d / "trait_events.tsv").write_text(_events_tsv(self.events))
    if "tree" in outputs:
        (d / "trait_tree.nwk").write_text(_trait_newick(self.complete_tree, self.node_values) + "\n")

zombi2.traits.Change dataclass

Change(time: float, kind: str, lineage: int, from_state: object, to_state: object)

A realized trait change — one entry of the event log, the trait twin of the genome level's :class:~zombi2.genomes.Event. On lineage lineage at time (crown-forward, the species-tree clock) the state went from from_state to to_state. kind is "on_branch" — a switch along a branch (an Mk transition) — "on_speciation" — a jump at a speciation node (from at_speciation; for a continuous trait from_state / to_state are the pre- and post-jump values) — or "root", one synthetic entry at the crown giving the initial state the run started in (from_state None, time the root's birth_time). That row is what lets the log stand on its own: the tree plus the root state plus the switches determines the trait on every lineage at every instant, so no separate driver file is needed.

Conditioning and joining

zombi2.joint.simulate_joint

simulate_joint(*, birth, death=0.0, trait=None, genome=None, n_extant=None, total_time=None, seed=None) -> JointResult

Grow a tree and the driver that drives its speciation, in one run (SPEC §2–4).

birth and death are rate specs (per lineage). Make either read the driver with mod.DrivenBy(source, mapping) — a live level name (not a filename) is what makes this joint rather than conditioned. Give exactly one driver:

• trait = traits.discrete(...) — a discrete trait drives speciation (BiSSE / MuSSE), read as mod.DrivenBy("trait", {"small": 1.0, "large": 2.0}). Driving both birth and death gives state-dependent λ and μ.
• genome = genomes.family(...) — gene content drives speciation (P(Species, Genomes)), read as the total gene count mod.DrivenBy("genomes:count", curve) or the presence of a named family mod.DrivenBy("genomes:toxin", {"present": 2.0, "absent": 1.0}) (declare it with family_names=["toxin"]).

::

joint.simulate_joint(
    birth  = 1.0 * mod.DrivenBy("genomes:toxin", {"present": 3.0, "absent": 1.0}),
    genome = genomes.family(origination=0.2, loss=0.1, family_names=["toxin"]),
    n_extant = 100, seed = 1)

The driver is an unexecuted process spec, grown with the tree. Stop at exactly n_extant living lineages (conditioned on survival — a birth-death tree can die out, so it restarts, advancing the same generator) or at total_time — give exactly one. Returns a :class:JointResult carrying the grown tree and the driver level (.trait or .genome). Deterministic given seed. Continuous trait→speciation (QuaSSE), clade drift (FromParent) combined with driving, and gene transfer in a joint run are not available.

Source code in zombi2/joint/__init__.py

def simulate_joint(*, birth, death=0.0, trait=None, genome=None, n_extant=None, total_time=None,
                   seed=None) -> JointResult:
    """Grow a tree **and** the driver that drives its speciation, in one run (SPEC §2–4).

    ``birth`` and ``death`` are rate specs (per lineage). Make either read the driver with
    ``mod.DrivenBy(source, mapping)`` — a **live level name** (not a filename) is what makes this
    *joint* rather than conditioned. Give **exactly one** driver:

    - ``trait = traits.discrete(...)`` — a discrete trait drives speciation (BiSSE / MuSSE), read as
      ``mod.DrivenBy("trait", {"small": 1.0, "large": 2.0})``. Driving both birth and death gives
      state-dependent λ *and* μ.
    - ``genome = genomes.family(...)`` — **gene content** drives speciation (``P(Species, Genomes)``),
      read as the total gene count ``mod.DrivenBy("genomes:count", curve)`` or the presence of a named
      family ``mod.DrivenBy("genomes:toxin", {"present": 2.0, "absent": 1.0})`` (declare it with
      ``family_names=["toxin"]``).

    ::

        joint.simulate_joint(
            birth  = 1.0 * mod.DrivenBy("genomes:toxin", {"present": 3.0, "absent": 1.0}),
            genome = genomes.family(origination=0.2, loss=0.1, family_names=["toxin"]),
            n_extant = 100, seed = 1)

    The driver is an **unexecuted** process spec, grown with the tree. Stop at exactly ``n_extant``
    living lineages (conditioned on survival — a birth-death tree can die out, so it restarts,
    advancing the same generator) **or** at ``total_time`` — give exactly one. Returns a
    :class:`JointResult` carrying the grown tree and the driver level (``.trait`` or ``.genome``).
    Deterministic given ``seed``. Continuous trait→speciation (QuaSSE), clade drift (``FromParent``)
    combined with driving, and gene transfer in a joint run are not available.
    """
    birth_rate = as_rate(birth, default_scope=PerLineage)
    death_rate = as_rate(death, default_scope=PerLineage)
    if (trait is None) == (genome is None):
        raise TypeError(
            "give exactly one driver: trait=traits.discrete(...) OR genome=genomes.family(...)."
        )
    # collect the DrivenBy sources on birth/death (a joint model's diversification must be per lineage)
    sources: list[str] = []
    for label, rate in (("birth", birth_rate), ("death", death_rate)):
        if not isinstance(rate.scope, PerLineage):
            raise ValueError(
                f"{label} has a {type(rate.scope).__name__} scope, but a joint diversification rate is "
                f"per lineage — drop the scope wrapper (per lineage is the default)."
            )
        for m in rate.modifiers:
            if isinstance(m, FromParent):
                raise ValueError(
                    f"{label} carries FromParent (clade drift); drift and a driven rate are not available "
                    f"together — use one or the other."
                )
            if isinstance(m, DrivenBy):
                if not isinstance(m.source, str):
                    raise TypeError(
                        f"{label} is driven by a {type(m.source).__name__} object, but a joint model "
                        f"drives from a live level *name* (a string, e.g. \"trait\" / \"genomes:count\"). "
                        f"A grown result object is conditioning — pass it to the target level's run."
                    )
                sources.append(m.source)
    if not sources:
        raise ValueError(
            "a joint model needs the driver to drive something: give birth (or death) a "
            "mod.DrivenBy(...). With neither driven, grow the two levels as independent runs instead."
        )
    # the driver spec must match the sources
    if trait is not None:
        if not isinstance(trait, DiscreteTrait):
            raise TypeError(
                "trait= must be traits.discrete(states=[...], switch=...) — a discrete process spec. "
                "Continuous trait→speciation (QuaSSE) is not available: a continuously varying "
                "rate needs thinning, which this exact race does not do."
            )
        bad = sorted({s for s in sources if s != "trait"})
        if bad:
            raise ValueError(
                f'with trait=, drive from the live trait — mod.DrivenBy("trait", ...); got source(s) '
                f"{bad}. (A filename source is conditioning, not a joint run.)"
            )
    else:
        if not isinstance(genome, FamilyGenome):
            raise TypeError("genome= must be genomes.family(...) — a family-genome process spec.")
        for s in sources:
            if s == _GENOME_COUNT:
                continue
            if s.startswith("genomes:"):
                name = s.split(":", 1)[1]
                if name not in genome.family_names:
                    raise ValueError(
                        f'DrivenBy("{s}", ...) names family {name!r}, but genomes.family was not '
                        f"declared with it — add family_names=[…, {name!r}]."
                    )
                continue
            raise ValueError(
                f'with genome=, drive from gene content — "genomes:count" or "genomes:<family>"; '
                f"got {s!r}."
            )
    if (n_extant is None) == (total_time is None):
        raise ValueError("give exactly one of n_extant or total_time")
    if n_extant is not None and (isinstance(n_extant, bool) or not isinstance(n_extant, int) or n_extant < 1):
        raise ValueError(f"n_extant must be a positive integer, got {n_extant!r}")
    if total_time is not None and (not isinstance(total_time, (int, float))
                                   or not math.isfinite(total_time) or total_time <= 0):
        raise ValueError(f"total_time must be a positive finite number, got {total_time!r}")

    rng = np.random.default_rng(seed)
    unique_sources = sorted(set(sources))

    def grow_once(target_n, tt) -> tuple[Tree, JointResult]:
        if trait is not None:
            tree, se, nv, te = _grow_joint(rng, birth_rate, death_rate, trait, target_n, tt)
            te.sort(key=lambda c: c.time)
            result = JointResult(SpeciesResult(tree, se, seed, []), seed,
                                 trait=TraitsResult(tree, nv, te, seed, kind="discrete"))
        else:
            tree, se, go, ge, fn = _grow_joint_genome(
                rng, birth_rate, death_rate, genome, unique_sources, target_n, tt)
            result = JointResult(SpeciesResult(tree, se, seed, []), seed,
                                 genome=FamilyGenomesResult(tree, go, ge, seed, fn))
        return tree, result

    if total_time is not None:
        return grow_once(None, total_time)[1]

    for _ in range(_MAX_ATTEMPTS):
        tree, result = grow_once(n_extant, None)
        if sum(1 for nd in tree.nodes.values() if nd.fate == "extant") == n_extant:
            return result
    raise RuntimeError(
        f"could not grow a tree to {n_extant} extant lineages in {_MAX_ATTEMPTS} attempts; "
        "birth must comfortably exceed death for large n_extant"
    )

zombi2.joint.JointResult dataclass

JointResult(species: SpeciesResult, seed: int | None, trait: TraitsResult | None = None, genome: FamilyGenomesResult | None = None)

What :func:simulate_joint returns — both grown levels of a joint run. species is the grown tree (a :class:~zombi2.species.SpeciesResult: complete_tree, extant_tree, the speciation/extinction events); the driver level that grew with it is either trait (a :class:~zombi2.traits.TraitsResult, for a trait→speciation run) or genome (a :class:~zombi2.genomes.FamilyGenomesResult, for a gene-content→speciation run) — exactly one is set. The tree is an output, grown by the driver it carries, so the levels share one complete_tree.

events property

events: list

The species events (speciation / extinction). The driver level's own events are trait.events / genome.events.

write

write(directory) -> None

Write both levels to directory: the species files (species_complete.nwk / species_extant.nwk / species_events.tsv) and the driver level's — for a trait, trait_values.tsv / trait_events.tsv / trait_tree.nwk; for a genome, genome_events.tsv / profiles.tsv.

Source code in zombi2/joint/__init__.py

def write(self, directory) -> None:
    """Write both levels to ``directory``: the species files (``species_complete.nwk`` /
    ``species_extant.nwk`` / ``species_events.tsv``) and the driver level's — for a trait,
    ``trait_values.tsv`` / ``trait_events.tsv`` / ``trait_tree.nwk``; for a genome,
    ``genome_events.tsv`` / ``profiles.tsv``."""
    self.species.write(directory, outputs=("complete", "extant", "events"))
    if self.trait is not None:
        self.trait.write(directory, outputs=("values", "events", "tree"))
    if self.genome is not None:
        self.genome.write(directory, outputs=("events", "profiles"))

Rates

Rates are written as expressions over scopes and modifiers; the same notation is shared by the Python API, the CLI and a --params file.

zombi2.rates.scope

Scope wrappers — the "per what?" of a rate (SPEC §5).

Every rate in ZOMBI2 is scope(base) × modifiers. A scope wrapper answers per what?: it tags a base rate with the unit it applies to, so the engine can turn that per-unit base into a total rate by multiplying by however many of the unit are present right now::

birth = scope.PerLineage(1.0)   # each lineage speciates at 1.0 -> total = 1.0 × (lineages alive)
birth = scope.Global(1.0)       # one shared budget for the whole tree -> total = 1.0 (constant)
loss  = scope.PerCopy(0.25)     # each gene copy is lost at 0.25 -> total = 0.25 × (copies present)

Wrappers wrap a base; they do not multiply — multiplying is what modifiers (zombi2.modifiers) do. The word "per" is reserved for these wrappers; a modifier never starts with "per".

There is deliberately no PerGenome: one genome lives in one lineage, so "per genome" is PerLineage.

A bare number (birth = 1.0) is coerced by each level to its natural default wrapper — species birth/death and gene origination default to per lineage, duplication/transfer/loss to per copy, substitution to per site. The wrappers here are the explicit override.

Scope dataclass

Scope(base: float)

A base rate tagged with the unit it applies to.

Abstract: use one of :class:Global, :class:PerLineage, :class:PerCopy, :class:PerSite, :class:PerChromosome.

total

total(**counts: float) -> float

The total rate given the current counts.

counts supplies the units in scope right now (lineages, copies, sites, chromosomes); each wrapper reads only the one it needs and ignores the rest. :class:Global reads none.

Source code in zombi2/rates/scope.py

def total(self, **counts: float) -> float:
    """The total rate given the current counts.

    ``counts`` supplies the units in scope right now (``lineages``, ``copies``,
    ``sites``, ``chromosomes``); each wrapper reads only the one it needs and
    ignores the rest. :class:`Global` reads none.
    """
    if self.unit is None:
        return self.base
    try:
        return self.base * counts[self.unit]
    except KeyError:
        raise KeyError(
            f"{type(self).__name__} needs a {self.unit!r} count; got {sorted(counts)}"
        ) from None

Global dataclass

Global(base: float)

Bases: Scope

One shared budget for the whole system: the total does not scale with anything.

Global (capitalised — global is a Python keyword) makes a process run at a constant total rate: linear growth, not exponential.

PerLineage dataclass

PerLineage(base: float)

Bases: Scope

Per lineage — the total scales with the number of lineages present.

The default for species birth/death and gene origination. Within a single genome there is one lineage, so this reads as a constant per-genome budget; across the species tree it is base × (lineages alive) (exponential diversification).

PerCopy dataclass

PerCopy(base: float)

Bases: Scope

Per gene copy — the total scales with family/genome size (duplication, transfer, loss).

A large family therefore turns over faster: base × (copies present).

PerSite dataclass

PerSite(base: float)

Bases: Scope

Per sequence site — the total scales with the number of sites (substitutions).

PerChromosome dataclass

PerChromosome(base: float)

Bases: Scope

Per chromosome — the total scales with the number of chromosomes (fission/fusion/loss).

zombi2.rates.modifiers

Rate modifiers — the context multipliers of a rate (SPEC §5).

Every rate in ZOMBI2 is scope(base) × modifiers. A modifier multiplies a rate by a dimensionless factor that depends on context — the current time, the standing diversity, the branch, the family, a driver level. Modifiers multiply (that is the whole difference from scope wrappers, which wrap), and the word "per" is reserved for scope, so a modifier never starts with "per".

You reach them through mod::

birth = 1.0 * mod.OnTime({0: 1.0, 3: 0.3})   # a skyline: 1.0, then 0.3 from time 3 on
birth = 1.0 * mod.OnTotalDiversity(cap=100)       # slows to 0 as diversity approaches 100

The deterministic modifiers (OnTime, OnTotalDiversity) compute their factor as a pure function of the context. The stochastic ones also carry a draw the engine drives with a random generator: FromParent (the rate drifts parent→child, via initial/descend), ByLineage (one i.i.d. draw per lineage) and ByFamily (one per family), the last two via draw.

Composition (*) belongs to the Rate module; a modifier here knows only how to produce its own factor, or its own draw.

Modifier

Base for rate modifiers.

A modifier reads the context keys it cares about (time, diversity, branch, family, …) and returns a dimensionless, non-negative multiplier; it ignores the rest. Abstract — use a subclass.

next_change

next_change(time: float) -> float

The next time strictly after time at which this modifier's factor changes on its own — a skyline breakpoint. inf if it never changes with time (the default; most modifiers change only at events, not autonomously).

Source code in zombi2/rates/modifiers.py

def next_change(self, time: float) -> float:
    """The next time strictly after ``time`` at which this modifier's factor changes on
    its own — a skyline breakpoint. ``inf`` if it never changes with time (the default;
    most modifiers change only at events, not autonomously)."""
    return math.inf

OnTime

OnTime(schedule: Mapping[float, float])

Bases: Modifier

The rate changes in time — a skyline / episodic schedule.

schedule maps each interval's start time to a relative factor::

OnTime({0: 1.0, 3: 0.3})   # factor 1.0 on [0, 3), then 0.3 from time 3 on

Factors are relative (dimensionless): on a base of 2.0 the schedule scales it. Before the earliest breakpoint the earliest factor applies (define the schedule from time 0 to avoid surprise).

Source code in zombi2/rates/modifiers.py

def __init__(self, schedule: Mapping[float, float]) -> None:
    steps = tuple(sorted((float(t), float(f)) for t, f in schedule.items()))
    if not steps:
        raise ValueError("OnTime needs a non-empty schedule, e.g. OnTime({0: 1.0, 3: 0.3})")
    for t, f in steps:
        if not math.isfinite(t):
            raise ValueError(f"OnTime schedule times must be finite, got {t!r}")
        if not math.isfinite(f) or f < 0:
            raise ValueError(f"OnTime factors must be finite and non-negative, got {f!r}")
    self._steps = steps

OnTotalDiversity dataclass

OnTotalDiversity(cap: float)

Bases: Modifier

The rate slows as standing diversity grows — diversity-dependence.

The factor falls linearly from 1 toward 0 as diversity rises to cap (a carrying capacity), and stays 0 beyond it: OnTotalDiversity(cap=100) halves the rate at 50 lineages and stops it at 100.

FromParent dataclass

FromParent(spread: float)

Bases: Modifier

The rate drifts along the tree — each lineage inherits its parent's rate times a random factor drawn at the split (geometric Brownian motion on the rate: clade drift at the species level, the autocorrelated clock at the sequence level). spread (σ) sets the drift width.

The per-split factor is lognormal, mean-corrected so E[factor] = 1. Without the correction the rate inflates down the tree (E[rate] ≈ e^{σ²/2} instead of 1) — a real historical bug. The draw logic (:meth:initial / :meth:descend) is driven by the engine, which threads each lineage's current factor and passes it back to :meth:factor as inherited.

initial

initial() -> float

The root's factor: 1.0 — the rate starts at its base.

Source code in zombi2/rates/modifiers.py

def initial(self) -> float:
    """The root's factor: 1.0 — the rate starts at its base."""
    return 1.0

descend

descend(parent_value: float, rng) -> float

A daughter's factor: the parent's, times one mean-corrected lognormal step.

Source code in zombi2/rates/modifiers.py

def descend(self, parent_value: float, rng) -> float:
    """A daughter's factor: the parent's, times one mean-corrected lognormal step."""
    sigma = self.spread
    return parent_value * math.exp(rng.normal(-0.5 * sigma * sigma, sigma))

factor

factor(*, inherited: float = 1.0, **_: float) -> float

The lineage's current factor — the engine threads it and passes it back as inherited.

Source code in zombi2/rates/modifiers.py

def factor(self, *, inherited: float = 1.0, **_: float) -> float:
    """The lineage's current factor — the engine threads it and passes it back as ``inherited``."""
    return inherited

ByLineage dataclass

ByLineage(spread: float, dist: str = 'lognormal')

Bases: Modifier

The rate varies independently from lineage to lineage — an uncorrelated ("relaxed") clock.

Each lineage draws one i.i.d. multiplier with no memory of its parent (contrast :class:FromParent, whose rate drifts parent→child). The draw is mean-corrected so E[factor] = 1 — without it the mean rate inflates down the tree (the historical lognormal-clock bug). spread (σ) sets the width; dist is "lognormal" (default; σ = the log-scale) or "gamma" (σ = the coefficient of variation) — the two agree to first order in σ.

At the sequence level this is the lineage clock: the engine draws one value per species lineage (via :meth:draw) and shares it across every gene family passing through that lineage. It is the lineage-twin of the genome level's ByFamily.

draw

draw(rng) -> float

One independent, mean-1 multiplier for a lineage. spread = 0 gives 1.0 (a strict clock).

Source code in zombi2/rates/modifiers.py

def draw(self, rng) -> float:
    """One independent, mean-1 multiplier for a lineage. ``spread = 0`` gives 1.0 (a strict clock)."""
    s = self.spread
    if s == 0.0:
        return 1.0
    if self.dist == "lognormal":
        return math.exp(rng.normal(-0.5 * s * s, s))     # mean-corrected lognormal
    return float(rng.gamma(1.0 / (s * s), s * s))        # mean-1 gamma, coefficient of variation = s

factor

factor(*, bylineage: float = 1.0, **_: float) -> float

The lineage's drawn factor — the engine threads it and passes it back as bylineage.

Source code in zombi2/rates/modifiers.py

def factor(self, *, bylineage: float = 1.0, **_: float) -> float:
    """The lineage's drawn factor — the engine threads it and passes it back as ``bylineage``."""
    return bylineage

ByFamily dataclass

ByFamily(spread: float, dist: str = 'lognormal')

Bases: Modifier

The rate varies independently from gene family to gene family.

The family-twin of :class:ByLineage, and the same i.i.d.-heterogeneity idea: each family draws one multiplier with no memory, mean-corrected so E[factor] = 1 — so widening spread spreads the families out without moving the average one off the base rate. dist is "lognormal" (default; σ = the log-scale) or "gamma" (σ = the coefficient of variation).

Where you put it decides what varies together. On a single rate, that rate varies by family on its own::

loss = 0.25 * mod.ByFamily(spread=0.5)      # a family that loses fast is not thereby
duplication = 0.2 * mod.ByFamily(spread=0.5)   # duplicating fast — independent draws

In the family-wide family_speed= slot, one draw scales every rate that family has, so a fast family is fast at everything::

simulate_genomes_family(tree, duplication=0.2, loss=0.25,
                        family_speed=mod.ByFamily(spread=0.5))

The two compose: a family-wide tempo, plus extra variation on one rate.

Not accepted on origination, which is the rate at which families are created — at the moment it is read there is no family to have drawn a factor for. The engine rejects it rather than quietly ignoring it.

draw

draw(rng) -> float

One independent, mean-1 multiplier for a family. spread = 0 gives 1.0 (no variation).

Source code in zombi2/rates/modifiers.py

def draw(self, rng) -> float:
    """One independent, mean-1 multiplier for a family. ``spread = 0`` gives 1.0 (no variation)."""
    s = self.spread
    if s == 0.0:
        return 1.0
    if self.dist == "lognormal":
        return math.exp(rng.normal(-0.5 * s * s, s))     # mean-corrected lognormal
    return float(rng.gamma(1.0 / (s * s), s * s))        # mean-1 gamma, coefficient of variation = s

factor

factor(*, byfamily: float = 1.0, **_: float) -> float

The family's drawn factor — the engine threads it and passes it back as byfamily.

Source code in zombi2/rates/modifiers.py

def factor(self, *, byfamily: float = 1.0, **_: float) -> float:
    """The family's drawn factor — the engine threads it and passes it back as ``byfamily``."""
    return byfamily

DrivenBy

DrivenBy(source: object, mapping: object)

Bases: Modifier

The rate is driven by another level — the one coupling mechanism (SPEC §2).

A coupling is Ch2's definition made literal: a parameter that reads its value from another level instead of a number you type. DrivenBy reads the driver's value on each lineage and multiplies the base rate by the mapped factor::

loss = 0.25 * mod.DrivenBy("habitat.tsv", {"aquatic": 3.0, "terrestrial": 1.0})
birth = 1.0 * mod.DrivenBy("trait", {"small": 1.0, "large": 2.0})   # a joint model

source says where the driver comes from, and that single choice splits conditioned from joint — the chapter's spine, can the driver be grown first?:

• a filename ("habitat.tsv") or a grown driver result (a discrete TraitsResult) — the driver was grown first and handed over (conditioned): two ordinary runs. The result object is the file's in-memory shortcut — same conditioning, no write/read step;
• a level name ("trait", "genomes:count") — the driver co-evolves in one run (joint): neither level can be grown first.

mapping says how the driver's value becomes the factor — a :class:~zombi2.rates.mapping.Table (a dict, for a discrete driver), a :class:~zombi2.rates.mapping.Curve (a callable, continuous), or a :class:~zombi2.rates.mapping.Scalar (a log-link coefficient); a raw dict / callable / number is coerced (:func:~zombi2.rates.mapping.as_mapping).

Like :class:FromParent (inherited) and :class:ByLineage (bylineage), DrivenBy reads a value the engine threads per lineage — here a drivers mapping {source: value} — and is otherwise dumb: it just maps the value to a factor. The engine owns where the value comes from (a file it loaded, or the live level growing beside the tree) and when it changes (a discrete driver switches mid-branch, so the engine steps its Gillespie at each switch); a rate reaching an engine that has not threaded its source gets a factor of 1.0 (inert). DrivenBy targets a rate (a "how often") and multiplies; it does not drive a value, such as an OU optimum.

Source code in zombi2/rates/modifiers.py

def __init__(self, source: object, mapping: object) -> None:
    from .mapping import as_mapping

    if isinstance(source, str):
        if not source.strip():
            raise ValueError("DrivenBy source must be a non-empty string (a filename or level name)")
        self.key: object = source                # a string source is its own context key
    else:
        self.key = id(source)                    # an in-memory driver result (conditioning): key by identity
    self.source = source
    self.mapping = as_mapping(mapping)

factor

factor(*, drivers: Mapping | None = None, **_: float) -> float

The mapped multiplier for this lineage's driver value — the engine threads the value under drivers[key] (key is the source string, or the identity of an in-memory driver). No drivers (or this source absent) ⇒ 1.0, so an unthreaded rate is inert (the engine is responsible for supplying the value where the coupling is supported).

Source code in zombi2/rates/modifiers.py

def factor(self, *, drivers: Mapping | None = None, **_: float) -> float:
    """The mapped multiplier for this lineage's driver value — the engine threads the value under
    ``drivers[key]`` (``key`` is the source string, or the identity of an in-memory driver). No
    ``drivers`` (or this source absent) ⇒ 1.0, so an unthreaded rate is inert (the engine is
    responsible for supplying the value where the coupling is supported)."""
    if drivers is None:
        return 1.0
    value = drivers.get(self.key)
    if value is None:
        return 1.0
    return self.mapping.multiplier(value)