zombi2.species

Level 1: the species tree every other level lives on. One forward engine, shaped by scopes and modifiers rather than by a model zoo.

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, or a scope with verbs chained onto it); the default scope is per lineage (each lineage speciates/dies at the base rate, so the tree grows exponentially). Yule = death=0.

Rates that vary from lineage to lineage. PerLineage(1.0).varying_among('lineages', Drift(LogNormal(0.0, σ))) is inherited variation — a daughter starts from its parent's rate and is nudged at the split, so fast clades stay fast (clade drift; the literature's ClaDS). PerLineage(1.0).varying_among('lineages', LogNormal(0.0, σ)) is independent variation — every lineage draws its own multiplier with no memory of its parent (relaxed rates), which is the null to compare drift against: the same amount of rate heterogeneity, none of it heritable, so the tree-shape signature that heritability leaves (lopsidedness — fast clades hoarding the tips) is absent. Both draws are mean-corrected, so widening the law's σ spreads lineages out without moving the average one off the base rate; both make the lineage that speciates or dies drawn weighted by its own rate; and both must be counted per lineage. They answer the same question and a rate carrying both is refused.

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. Deterministic given seed.

Where n_extant puts the present. The run stops the first moment n_extant lineages are alive together, then draws one more waiting time and places the present where that next event would have fired, without applying it — so the two newest tips get a real branch length rather than a zero-length one. Under pure birth this is exactly the general sampling approach (Hartmann, Wong & Stadler 2010): the tree is the process observed at an instant drawn uniformly over the time it spends holding n lineages. With extinction it is not. It is a first hitting rule — the run stops the first time it touches n, so an interval at n lineages reached by falling back from n+1 is never sampled — and the trees are correspondingly shallower than the birth-death process conditioned on n tips. The gap grows with turnover and shrinks with n: it is within noise at death=0, around a tenth of the tree height at n=10 with death/birth=0.4, and roughly a third to a half at n=10 with death/birth=0.8; by n=50 at moderate turnover it is back in the noise. If you publish trees grown this way, say which rule made them — a rate estimator applied to them will otherwise look broken for reasons that are not its fault.

total_time is not conditioned on survival: it can die out, and then it raises rather than handing back a tree with no present. Looping over seeds and skipping the failures is survival conditioning by another name, and changes the distribution of everything downstream.

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, or a scope with verbs chained onto it); the
    default scope is **per lineage** (each lineage speciates/dies at
    the base rate, so the tree grows exponentially). Yule = ``death=0``.

    **Rates that vary from lineage to lineage.**
    ``PerLineage(1.0).varying_among('lineages', Drift(LogNormal(0.0, σ)))`` is *inherited*
    variation — a daughter starts from its parent's rate and is nudged at the split, so fast clades
    stay fast (clade drift; the literature's ClaDS).
    ``PerLineage(1.0).varying_among('lineages', LogNormal(0.0, σ))`` is
    *independent* variation — every lineage draws its own multiplier with no memory of its parent
    (*relaxed* rates), which is the null to compare drift against: the same amount of rate
    heterogeneity, none of it heritable, so the tree-shape signature that heritability leaves
    (lopsidedness — fast clades hoarding the tips) is absent. Both draws are mean-corrected, so
    widening the law's σ spreads lineages out without moving the average one off the base rate; both
    make the lineage that speciates or dies drawn weighted by its own rate; and both must be counted
    per lineage. They answer the same question and a rate carrying both is refused.

    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``. Deterministic given
    ``seed``.

    **Where ``n_extant`` puts the present.** The run stops the first moment ``n_extant`` lineages
    are alive together, then draws one more waiting time and places the present where that next
    event *would* have fired, without applying it — so the two newest tips get a real branch length
    rather than a zero-length one. Under pure birth this is exactly the general sampling approach
    (Hartmann, Wong & Stadler 2010): the tree is the process observed at an instant drawn uniformly
    over the time it spends holding n lineages. **With extinction it is not.** It is a *first
    hitting* rule — the run stops the first time it touches n, so an interval at n lineages reached
    by falling back from n+1 is never sampled — and the trees are correspondingly shallower than the
    birth-death process conditioned on n tips. The gap grows with turnover and shrinks with n: it is
    within noise at ``death=0``, around a tenth of the tree height at n=10 with ``death/birth=0.4``,
    and roughly a third to a half at n=10 with ``death/birth=0.8``; by n=50 at moderate turnover it
    is back in the noise. If you publish trees grown this way, say which rule made them — a rate
    estimator applied to them will otherwise look broken for reasons that are not its fault.

    ``total_time`` is not conditioned on survival: it can die out, and then it **raises** rather than
    handing back a tree with no present. Looping over seeds and skipping the failures is survival
    conditioning by another name, and changes the distribution of everything downstream.

    ``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)):
        scope = rate.scope
        assert scope is not None            # `as_rate` filled the level's default just above
        # 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 issubclass(scope, IMPLEMENTED_SCOPES):
            raise ValueError(
                f"{label} has a {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 m.reads == (DRAWN, "families"):
                # not a missing feature: there is nothing here for it to mean
                raise ValueError(
                    f"{label} varies among families, but a species tree has no gene families — "
                    f"varying_among('families', ...) belongs on a genomes rate. For per-lineage "
                    f"heterogeneity here use varying_among('lineages', LogNormal(0.0, 0.5)) "
                    f"(independent) or varying_among('lineages', Drift(LogNormal(0.0, 0.5))) "
                    f"(inherited)."
                )
            if not is_implemented(m, IMPLEMENTED_MODIFIERS, "species"):
                # A driven rate is missing from that list for a reason worth saying, and only a
                # driven one: a driver of speciation cannot be grown before the tree, because the
                # tree is what it would have to grow on, so the model is a joint run rather than a
                # species run (SPEC §2–4). Said on any other rejection it would only mislead.
                driven = m.reads is not None and m.reads[0] == DRIVEN
                raise ValueError(
                    f"{label} carries {describe(m)}, which the species engine does not "
                    f"support. It takes changing_at (skyline), scaled_by(TotalDiversity(cap=...)) "
                    f"(diversity-dependent), varying_among('lineages', Drift(...)) (inherited rate "
                    f"drift, ClaDS) and varying_among('lineages', dist) (independent per-lineage "
                    f"rates)."
                    + (" A trait or gene content that drives speciation has to grow with the tree, "
                       "since the tree is what it would grow on, so the model is a joint run: "
                       "joint.simulate(species.birth_death(birth=PerLineage(1.0).scaled_by('trait', {...})), ...)."
                       if driven else "")
                )
        # SPEC §5: one memory structure per axis. A bare distribution has no memory and a Drift has a
        # continuous one, so a rate carrying both asks for a lineage's factor to be independent of its
        # parent's and inherited from it at once — there is no model there to implement, so say so
        # rather than silently letting whichever comes first win.
        check_one_memory(_per_lineage(rate), label=label, unit="lineages")
        if (per_lineage := _per_lineage(rate)) and scope is not PerLineage:
            raise ValueError(
                f"{label} carries {describe(per_lineage[0])} (a factor per lineage) but its scope "
                f"is {scope.__name__}; a rate that varies by lineage must be counted per "
                f"lineage — write PerLineage(...), or a bare number, which is per lineage here."
            )
    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, seed = stream("species", seed)     # own stream, and a drawn seed if none was given

    def _finish(tree: Tree, events: list[Event],
                rates: dict[int, tuple[float, float, float, int]]) -> SpeciesResult:
        # observe (sampling relabels survivors) then recover fossils along the grown branches
        alive = sum(1 for nd in tree.nodes.values() if nd.fate == "extant")
        _apply_sampling(tree, sampling, rng)
        # Sampling can take none of them, and then the run has no present — the same dead end the
        # extinction guard above refuses, reached by the other road. Refusing here too keeps the two
        # "nothing observed" outcomes consistent: neither hands back a result whose extant tree is
        # None for a downstream level to trip over. Only reachable with sampling < 1, since both
        # callers guarantee a survivor before this point.
        if not any(nd.fate == "extant" for nd in tree.nodes.values()):
            raise RuntimeError(
                f"sampling={sampling:g} observed none of the {alive} survivor"
                f"{'' if alive == 1 else 's'}, so the run has no present to grow a genome, sequence "
                f"or trait along. This is the sampling process, not a bad parameter — it has "
                f"probability {(1 - sampling) ** alive:.3g} here — so raise sampling, ask for more "
                f"survivors, or draw another seed.")
        return SpeciesResult(tree, events, seed, _recover_fossils(tree, fossils, rng), rates,
                             {"birth": birth_rate, "death": death_rate})

    if total_time is not None:
        tree, events, rates = _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, rates)

    for _ in range(_MAX_ATTEMPTS):
        tree, events, rates = _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, rates)
    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(), rates_at_birth: dict[int, tuple[float, float, float, int]] = dict(), _rates: dict = dict())

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.

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, which simulate_species_tree refuses to return: a run with no present raises there instead, so a result that came from it always has one.

lineage_rates

lineage_rates(kind: str = 'birth') -> dict[str, float]

{lineage: rate} — the rate that lineage itself ran under, in events per unit time.

A rate that varies among lineages gives every lineage a factor of its own, and the tree records only what happened, not the rate it happened under: two runs with the same shape can come from very different rates. This is that number, for the lineage rather than for the run — the base with this lineage's own factors applied — so a branch can be coloured by it.

Taken at the lineage's birth, which is when its factors are drawn. A rate that also depends on time or on standing diversity keeps moving along the branch afterwards; this is its value at the start of it. A rate that varies among lineages in no way gives the same number for every lineage, which is correct rather than useless: it is that rate.

Lineages are named as every output file names them (n12), so the result drops straight into a plot keyed by tip name.

Source code in zombi2/species/__init__.py

def lineage_rates(self, kind: str = "birth") -> dict[str, float]:
    """``{lineage: rate}`` — the rate that lineage itself ran under, in events per unit time.

    A rate that varies among lineages gives every lineage a factor of its own, and the tree
    records only what happened, not the rate it happened under: two runs with the same shape can
    come from very different rates. This is that number, for the lineage rather than for the run
    — the base with this lineage's own factors applied — so a branch can be coloured by it.

    Taken **at the lineage's birth**, which is when its factors are drawn. A rate that also
    depends on time or on standing diversity keeps moving along the branch afterwards; this is
    its value at the start of it. A rate that varies among lineages in no way gives the same
    number for every lineage, which is correct rather than useless: it is that rate.

    Lineages are named as every output file names them (``n12``), so the result drops straight
    into a plot keyed by tip name.
    """
    if kind not in ("birth", "death"):
        raise ValueError(f"kind must be 'birth' or 'death', got {kind!r}")
    rate = self._rates.get(kind)
    if rate is None:                       # a result built by hand, or unpickled from an old run
        raise ValueError(f"this result carries no {kind} rate to evaluate")
    label = self.complete_tree.labels()
    j = 0 if kind == "birth" else 1
    return {label[i]: rate.effective(carried_factor=r[j], time=r[2], diversity=r[3], lineages=1)
            for i, r in self.rates_at_birth.items() if i in label}

summary

summary() -> dict

What this run produced, as a plain dict — the payload of species_summary.json.

Counts, not parameters: the log already says what was asked for. The realised rates are here because they are the cheapest check anyone can make on a tree — events divided by the exposure that generated them, which is what a declared per-lineage rate means.

Source code in zombi2/species/__init__.py

def summary(self) -> dict:
    """What this run produced, as a plain dict — the payload of ``species_summary.json``.

    Counts, not parameters: the log already says what was asked for. The realised rates are here
    because they are the cheapest check anyone can make on a tree — events divided by the exposure
    that generated them, which is what a declared per-lineage rate means."""
    nodes = self.complete_tree.nodes
    tips = [n for n in nodes.values() if not n.children]
    extant = self.complete_tree.extant_leaves()
    speciations = sum(1 for e in self.events if e.kind == "speciation")
    extinctions = sum(1 for e in self.events if e.kind == "extinction")
    # total branch length: every node's own branch, which is the exposure a per-lineage rate ran on
    exposure = sum(n.end_time - n.birth_time for n in nodes.values())
    height = max(nodes[i].end_time for i in extant) if extant else None
    root = nodes[self.complete_tree.root]
    return {
        "level": "species",
        "seed": self.seed,
        "tips": {"extant": len(extant), "extinct": len(self.complete_tree.extinct_leaves()),
                 "unsampled": len(self.complete_tree.unsampled_leaves()), "total": len(tips)},
        "nodes": len(nodes),
        "events": {"speciation": speciations, "extinction": extinctions},
        "fossils": len(self.fossils),
        "tree": {"height": height, "stem_length": root.end_time - root.birth_time,
                 "total_branch_length": exposure},
        # events per lineage per unit time, as declared rates are counted. A sanity check, not a
        # parameter: it is what the run realised, which a conditioned stop condition can bias.
        "realised_rates": {
            "birth": round(speciations / exposure, 6) if exposure else None,
            "death": round(extinctions / exposure, 6) if exposure else None},
    }

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, time · kind · parents · children), "fossils" → species_fossils.tsv (if any recovered), "fates" → species_fates.tsv (each tip's resolved fate), "summary" → species_summary.json (the summary() payload).

species_events.tsv names the lineages an event consumed and the lineages it produced, the same parents / children pair every event file uses: a speciation row is one parent and its two children (;-packed), an extinction row is the dying lineage as the parent with no children.

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, ``time`` · ``kind`` · ``parents`` · ``children``),
    ``"fossils"`` → ``species_fossils.tsv`` (if any recovered),
    ``"fates"`` → ``species_fates.tsv`` (each tip's resolved fate),
    ``"summary"`` → ``species_summary.json`` (the ``summary()`` payload).

    ``species_events.tsv`` names the lineages an event consumed and the lineages it produced, the
    same ``parents`` / ``children`` pair every event file uses: a ``speciation`` row is one parent
    and its two children (``;``-packed), an ``extinction`` row is the dying lineage as the parent
    with no children.

    ``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", encoding="utf-8")
    if "extant" in outputs and self.extant_tree is not None:
        (d / "species_extant.nwk").write_text(self.extant_tree.to_newick() + "\n", encoding="utf-8")
    name = self.complete_tree.labels()
    if "events" in outputs:
        # parents / children, not lineage / children: the lineage an event consumed IS its parent
        # (the one that split, or the one that died), so one column pair reads right for both kinds
        rows = ["time\tkind\tparents\tchildren"]
        for e in self.events:
            kids = ";".join(name[c] for c in e.children) if e.children else ""
            rows.append(f"{e.time:.6g}\t{e.kind}\t{name[e.node]}\t{kids}")
        (d / "species_events.tsv").write_text("\n".join(rows) + "\n", encoding="utf-8")
    if "fossils" in outputs and self.fossils:
        # fossils are drawn along every branch, a surviving lineage's as readily as an extinct one's
        rows = ["lineage\ttime"] + [f"{name[i]}\t{t:.6g}" for i, t in self.fossils]
        (d / "species_fossils.tsv").write_text("\n".join(rows) + "\n", encoding="utf-8")
    if "summary" in outputs:
        write_summary(d / "species_summary.json", self.summary())
    if "fates" in outputs:
        # one row per tip (extant / extinct / unsampled); internal nodes are always speciations
        rows = ["lineage\tfate"]
        for i in sorted(self.complete_tree.leaves()):
            rows.append(f"{name[i]}\t{self.complete_tree.nodes[i].fate}")
        (d / "species_fates.tsv").write_text("\n".join(rows) + "\n", encoding="utf-8")

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.

Trees

The tree object itself, and the readers and shape helpers that work on it. These live in zombi2.tree — from zombi2.tree import read_newick — and are documented here because the species level is the one that grows a tree; every other level takes one. read_newick reads a ZOMBI2 tree or an external one, so a genome run can start from a published phylogeny.

zombi2.tree.Tree dataclass

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

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

labels

labels() -> dict[int, str]

{node id: its written name} — n<id>, or e<id> for a lineage that went extinct.

The tree is the only thing that knows a lineage's fate, so this is where a run's names come from: a writer builds the map once and every id it prints goes through it. Not cached, because fates are assigned after construction when a tree is read back (read_newick()), and a map frozen before that would name every tip n.

Source code in zombi2/tree.py

def labels(self) -> dict[int, str]:
    """``{node id: its written name}`` — ``n<id>``, or ``e<id>`` for a lineage that went extinct.

    The tree is the only thing that knows a lineage's fate, so this is where a run's names come
    from: a writer builds the map once and every id it prints goes through it. Not cached, because
    fates are assigned after construction when a tree is read back (`read_newick()`), and a map
    frozen before that would name every tip ``n``."""
    return {i: node_label(i, n.fate) for i, n in self.nodes.items()}

leaves

leaves() -> list[int]

Every lineage with no descendants — extant and extinct.

Source code in zombi2/tree.py

def leaves(self) -> list[int]:
    """Every lineage with no descendants — extant **and** extinct."""
    return [i for i, n in self.nodes.items() if not n.children]

extant_leaves

extant_leaves() -> list[int]

The lineages alive at the present. (A tip list, not a tree — the pruned survivors' tree is SpeciesResult.extant_tree.)

Source code in zombi2/tree.py

def extant_leaves(self) -> list[int]:
    """The lineages alive at the present. (A *tip* list, not a tree — the pruned survivors' tree
    is `SpeciesResult.extant_tree`.)"""
    return [i for i, n in self.nodes.items() if n.fate == "extant"]

extinct_leaves

extinct_leaves() -> list[int]

The lineages that died before the present.

Source code in zombi2/tree.py

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

unsampled_leaves

unsampled_leaves() -> list[int]

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_leaves(self) -> list[int]:
    """Survivors not observed under incomplete ``sampling`` — kept in the complete tree (told
    apart by their fate) but pruned from the extant tree."""
    return [i for i, n in self.nodes.items() if n.fate == "unsampled"]

to_newick

to_newick(*, precision: int | None = None) -> 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>, or e<id> for a lineage that went extinct (see node_label()).

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 read_newick() reads it back.

precision is the number of significant digits each branch length is written to. None (the default) writes the shortest string that reads back as exactly the same float, so a tree written and re-read is the tree you had.

That exactness is not cosmetic: the CLI hands a tree between levels through this file, so a rounded length is a different tree. At the old fixed 12 digits every branch shifted by about 2e-12 on the round trip, which moved every downstream Gillespie waiting time — and zombi2 genomes --seed 7 and simulate_genomes_family(sp, seed=7) then produced different histories from the same tree and the same seed. Both were valid draws, but a seed that means one run through Python and another through the CLI is not a seed anyone can publish. Writing lengths in full costs a few bytes a branch and makes the two front doors the same run.

Digits also matter downstream. Seven — the default before 12 — was not enough: a tip's depth is a sum of branch lengths, so rounding accumulates down the path, and on a 40-tip tree of height 4 two tips written at 7 digits came out about 1e-6 apart. That is far above the tolerance ape::is.ultrametric() allows (~1e-8), so an ultrametric tree — which every extant tree from a dated run is, to 1e-16 in memory — was rejected by the first thing anyone does with it in R. Pass precision= to go back to fixed significant digits for a smaller file; the tree it writes no longer round-trips exactly.

Source code in zombi2/tree.py

def to_newick(self, *, precision: int | None = None) -> 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>``, or ``e<id>`` for a lineage that went extinct (see `node_label()`).

    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 `read_newick()` reads it back.

    ``precision`` is the number of **significant digits** each branch length is written to.
    ``None`` (the default) writes the shortest string that reads back as *exactly* the same
    float, so a tree written and re-read is the tree you had.

    That exactness is not cosmetic: the CLI hands a tree between levels through this file, so a
    rounded length is a different tree. At the old fixed 12 digits every branch shifted by about
    2e-12 on the round trip, which moved every downstream Gillespie waiting time — and
    ``zombi2 genomes --seed 7`` and ``simulate_genomes_family(sp, seed=7)`` then produced
    *different* histories from the same tree and the same seed. Both were valid draws, but a seed
    that means one run through Python and another through the CLI is not a seed anyone can
    publish. Writing lengths in full costs a few bytes a branch and makes the two front doors the
    same run.

    Digits also matter downstream. Seven — the default before 12 — was not enough: a tip's
    *depth* is a sum of branch lengths, so rounding accumulates down the path, and on a 40-tip
    tree of height 4 two tips written at 7 digits came out about 1e-6 apart. That is far above the
    tolerance ``ape::is.ultrametric()`` allows (~1e-8), so an ultrametric tree — which every
    extant tree from a dated run is, to 1e-16 in memory — was rejected by the first thing anyone
    does with it in R. Pass ``precision=`` to go back to fixed significant digits for a smaller
    file; the tree it writes no longer round-trips exactly."""

    name = self.labels()
    # repr() is Python's shortest round-tripping float form: float(repr(x)) == x, always
    num = repr if precision is None else (lambda x: f"{x:.{precision}g}")

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

    root = self.nodes[self.root]
    stem = root.end_time - root.birth_time
    if not root.children:
        return f"{name[self.root]}:{num(stem)};"
    return f"({','.join(emit(c) for c in root.children)}){name[self.root]}:{num(stem)};"

zombi2.tree.Node dataclass

Node(id: int, parent: int | None, birth_time: float, end_time: float = math.inf, children: tuple[int, int] | tuple[()] = (), 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.

end_time defaults to inf — has not ended yet, which is the state a lineage is in while the engine is still growing it. Every node of a finished Tree has a real end, so the toolkit reads it as the number it is; a stray inf that escaped would surface as an infinite branch length rather than as a None propagating quietly through the arithmetic.

is_leaf property

is_leaf: bool

No descendants — a tip. (GeneNode spells this the same way; the two node types diverge in how they store children, but "is this a tip" reads identically on both.)

zombi2.tree.prune

prune(tree: Tree, keep: str = 'extant', *, tips: 'set[int] | None' = None) -> 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, is not built: simulate_species_tree reports fossils as (lineage, time) pairs rather than as taxa, so there are no sampled ancestors to keep.

tips keeps a named set of leaves instead, whatever their fate — which is the same operation on a different question: "the tree of the survivors" against "the tree of these taxa". Comparing a gene tree to the species tree needs the second, because a family present in part of the tree can only be judged against the part it occupies; without it zombi2 tools treedist could score nothing but a universal single-copy family. The branch lengths merge across the suppressed nodes exactly as they do for the extant tree, so the pruned tree is a real dated tree and a length-aware metric means something on it.

Source code in zombi2/tree.py

def prune(tree: Tree, keep: str = "extant", *, tips: "set[int] | None" = None) -> 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, is not built: `simulate_species_tree` reports fossils as
    ``(lineage, time)`` pairs rather than as taxa, so there are no sampled ancestors to keep.

    ``tips`` keeps a **named set of leaves** instead, whatever their fate — which is the same
    operation on a different question: "the tree of the survivors" against "the tree of these taxa".
    Comparing a gene tree to the species tree needs the second, because a family present in part of
    the tree can only be judged against the part it occupies; without it `zombi2 tools treedist`
    could score nothing but a universal single-copy family. The branch lengths merge across the
    suppressed nodes exactly as they do for the extant tree, so the pruned tree is a real dated tree
    and a length-aware metric means something on it."""
    if tips is None and keep != "extant":
        raise ValueError(
            f"keep must be 'extant', got {keep!r}: the sampled (fossil) tree is not built — fossils "
            f"are reported as (lineage, time) pairs, not as taxa on a tree.")
    nodes = tree.nodes
    #: whether a leaf is one of the ones being kept — the only thing the two modes differ in
    def _keep_leaf(i: int) -> bool:
        return nodes[i].fate == "extant" if tips is None else i in tips
    surviving: dict[int, bool] = {}
    for i in sorted(nodes, reverse=True):  # children have higher ids → processed before parents
        nd = nodes[i]
        surviving[i] = _keep_leaf(i) if not nd.children 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 not nd.children 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 (not nodes[i].children and _keep_leaf(i)) 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, (), nodes[i].fate)
        if p is None:
            ext_root = i

    def kept_children(i: int) -> tuple[int, int] | None:
        """``i``'s kept children, **in the order the complete tree had them**.

        Descends through the surviving unifurcations a pruned sibling leaves behind, taking the
        original children left to right. Rebuilding them by node id instead — which this did — is
        stable but not *consistent*: ids are assigned in birth order, and a suppressed child is
        replaced by a descendant whose id can be far larger than its sibling's, so the pair comes
        out swapped. The complete tree and its extant tree then draw the same clade on opposite
        sides, and any figure showing both, or any reader joining them by position, disagrees with
        itself."""
        out: list[int] = []
        stack: list[int] = list(reversed(nodes[i].children or ()))
        while stack:
            c = stack.pop()
            if not surviving[c]:
                continue
            if c in kept:
                out.append(c)
            else:                                   # a unifurcation: its kept descendants stand in
                stack.extend(reversed(nodes[c].children or ()))
        if not out:
            return None
        if len(out) != 2:                       # `kept` is exactly the leaves and the real splits
            raise AssertionError(f"node {i} kept {len(out)} children; a pruned tree is bifurcating")
        return out[0], out[1]

    for i in kept:
        children = kept_children(i)
        if children is not None:
            new[i].children = children

    assert ext_root is not None                 # `kept` is non-empty, so its root was found above
    return Tree(new, ext_root)

zombi2.tree.read_newick

read_newick(newick: str, *, tip_fates: dict[str, str] | None = None, assume_extant: bool = False) -> tuple[Tree, dict[int, str]]

Parse a Newick string into a complete Tree and a name-map {id: user label}.

This is how the CLI loads a species tree back for the downstream levels. Branch lengths are read as durations: the root sits at time 0 and each node's birth_time is its parent's end_time, so end_time - birth_time is the parsed length. Two kinds of tree are accepted, told apart by the labels:

• a ZOMBI complete tree (every node — internal ones too — is n<id>, or e<id> for a lineage that died, as to_newick writes it): the ids come from the labels, and the name-map is empty (the labels are the ids). Fate comes from tip_fates when given — the run's species_fates.tsv, keyed by the same label — which is authoritative; without it an e<id> label is the fate, and only a tree carrying no e labels at all falls back to depth (a leaf is "extinct" if it ends before the tree's greatest depth, else "extant"). Neither fallback can recover an "unsampled" tip.

• any external tree (leaves named freely, internal nodes usually unlabelled): fresh ids are minted in traversal order (root 0, parents before children), the original labels are returned as the name-map ({minted id: user label}), and fates depend on whether the tree is ultrametric (all root-to-tip depths equal, within 1e-6 × height):

• ultrametric → the tips are contemporaneous, so every tip is "extant" (observed);

• not ultrametric → the differing tip depths could mean extinct lineages or early samples, which ZOMBI cannot tell apart, so it refuses to guess: pass tip_fates — a {tip label: "extant" | "extinct" | "unsampled"} map covering every tip — or a ValueError is raised. (The CLI fills tip_fates from --tip-fates FILE, which reads the same format a species run writes to species_fates.tsv.)

A root branch length is read when present — to_newick writes one, so a ZOMBI tree round-trips with its stem intact. External trees usually have none, and then the root gets zero duration and the tree starts at its crown, which is all the file says.

The .nwk records only branch lengths, so "unsampled" fate cannot be read from it alone — pass tip_fates (the run's species_fates.tsv) to recover it; without one a survivor reads back "extant", which is still fine for evolving genomes/traits along the tree.

Only bifurcating trees are supported (an internal node with other than two children raises).

Source code in zombi2/tree.py

def read_newick(newick: str, *, tip_fates: dict[str, str] | None = None,
                assume_extant: bool = False) -> tuple[Tree, dict[int, str]]:
    """Parse a Newick string into a complete `Tree` and a name-map ``{id: user label}``.

    This is how the CLI loads a species tree back for the downstream levels. Branch lengths are read
    as **durations**: the root sits at time 0 and each node's ``birth_time`` is its parent's
    ``end_time``, so ``end_time - birth_time`` is the parsed length. Two kinds of tree are accepted,
    told apart by the labels:

    - a **ZOMBI complete tree** (every node — internal ones too — is ``n<id>``, or ``e<id>`` for a
      lineage that died, as ``to_newick`` writes it): the ids come from the labels, and the name-map
      is empty (the labels *are* the ids). Fate comes from ``tip_fates`` when given — the run's
      ``species_fates.tsv``, keyed by the same label — which is authoritative; without it an
      ``e<id>`` label *is* the fate, and only a tree carrying no ``e`` labels at all falls back to
      depth (a leaf is ``"extinct"`` if it ends before the tree's greatest depth, else ``"extant"``).
      Neither fallback can recover an ``"unsampled"`` tip.
    - any **external tree** (leaves named freely, internal nodes usually unlabelled): fresh ids are
      minted in traversal order (root 0, parents before children), the original labels are returned as
      the **name-map** (``{minted id: user label}``), and fates depend on whether the tree is
      **ultrametric** (all root-to-tip depths equal, within ``1e-6 × height``):

      - **ultrametric** → the tips are contemporaneous, so **every tip is ``"extant"``** (observed);
      - **not ultrametric** → the differing tip depths could mean extinct lineages *or* early
        samples, which ZOMBI cannot tell apart, so it **refuses to guess**: pass ``tip_fates`` — a
        ``{tip label: "extant" | "extinct" | "unsampled"}`` map covering every tip — or a
        `ValueError` is raised. (The CLI fills ``tip_fates`` from ``--tip-fates FILE``, which
        reads the same format a species run writes to ``species_fates.tsv``.)

    A root branch length is read when present — ``to_newick`` writes one, so a ZOMBI tree round-trips
    with its stem intact. External trees usually have none, and then the root gets zero duration and
    the tree starts at its crown, which is all the file says.

    The ``.nwk`` records only branch lengths, so ``"unsampled"`` fate cannot be read from it alone —
    pass ``tip_fates`` (the run's ``species_fates.tsv``) to recover it; without one a survivor reads
    back ``"extant"``, which is still fine for evolving genomes/traits along the tree.

    Only bifurcating trees are supported (an internal node with other than two children raises).
    """
    s = newick.strip().rstrip(";").strip()
    if not s:
        raise ValueError("empty Newick string — is the tree file empty?")
    i = 0

    def skip_ws() -> None:
        # whitespace (incl. the newlines of a line-wrapped file) is insignificant between tokens
        nonlocal i
        while i < len(s) and s[i].isspace():
            i += 1

    def read_name() -> str:
        # a quoted label ('...' / "...") is taken verbatim (a doubled quote unwrapped to one); an
        # unquoted label runs to the next whitespace or structural char, so stray whitespace never leaks
        nonlocal i
        if i < len(s) and s[i] in "'\"":
            quote = s[i]
            i += 1
            chars: list[str] = []
            while i < len(s):
                if s[i] == quote:
                    if i + 1 < len(s) and s[i + 1] == quote:
                        chars.append(quote)
                        i += 2
                        continue
                    i += 1
                    break
                chars.append(s[i])
                i += 1
            return "".join(chars)
        start = i
        while i < len(s) and s[i] not in ",():;" and not s[i].isspace():
            i += 1
        return s[start:i]

    # parse into a lightweight tree of (name, length, children) so we can decide ids/fates in a
    # second pass, once we know whether every node is ``n<id>``-labelled and the tree's depth
    class _P:
        __slots__ = ("name", "length", "children")

        def __init__(self, name, length, children):
            self.name, self.length, self.children = name, length, children

    def parse() -> _P:
        nonlocal i
        children: list[_P] = []
        skip_ws()
        if i < len(s) and s[i] == "(":
            i += 1
            while True:
                children.append(parse())
                skip_ws()
                if i >= len(s):
                    raise ValueError("malformed Newick: unbalanced parentheses — the string ended "
                                     "while a clade was still open (is the tree file truncated?)")
                if s[i] == ",":
                    i += 1
                elif s[i] == ")":
                    i += 1
                    break
                else:
                    raise ValueError(f"malformed Newick: expected ',' or ')' after a clade at "
                                     f"position {i}, got {s[i]!r}")
        skip_ws()
        name = read_name()
        length = 0.0
        skip_ws()
        if i < len(s) and s[i] == ":":
            i += 1
            skip_ws()
            start = i
            while i < len(s) and s[i] not in ",():;" and not s[i].isspace():
                i += 1
            try:
                length = float(s[start:i])
            except ValueError:
                raise ValueError(f"malformed Newick: expected a branch length after ':' at position "
                                 f"{start}, got {s[start:i]!r}") from None
        if children and len(children) != 2:
            raise ValueError(f"only bifurcating trees are supported: node {name or '(unnamed)'!r} has "
                             f"{len(children)} children (collapse polytomies / unifurcations first)")
        return _P(name, length, children)

    root_p = parse()
    skip_ws()
    if i < len(s):
        raise ValueError(f"malformed Newick: unexpected trailing text at position {i}: {s[i:]!r}")

    # ZOMBI complete tree ⟺ every node carries an ``n<id>`` label; then ids come from the labels,
    # otherwise we mint them ourselves and call every tip extant.
    all_labelled = True

    def _scan(p: _P) -> None:
        nonlocal all_labelled
        if not _ZOMBI_LABEL.match(p.name):
            all_labelled = False
        for c in p.children:
            _scan(c)

    _scan(root_p)

    # An external tree's tip labels are the join back to the caller's own taxa — the CLI writes them
    # to ``names.tsv`` — so two tips sharing a label make that table ambiguous (``n2→A`` beside
    # ``n4→A``) and every downstream merge by taxon name silently duplicates or drops rows. The
    # simulation itself would be fine, which is what makes this worth refusing rather than warning
    # about: nothing downstream would ever look wrong. A ZOMBI tree cannot reach this — its labels
    # are ids, and repeats are caught as duplicate ids below.
    if not all_labelled:
        seen: set[str] = set()
        repeated: list[str] = []
        stack = [root_p]
        while stack:
            p = stack.pop()
            if p.children:
                stack.extend(p.children)
            elif p.name:
                (repeated.append(p.name) if p.name in seen else seen.add(p.name))
        if repeated:
            raise ValueError(
                f"duplicate tip label(s) in the Newick: {', '.join(sorted(set(repeated)))} — tip "
                f"labels are how results join back to your taxa, so they must be unique. Rename the "
                f"repeats (a species appearing twice is usually an export or hand-edit slip).")

    nodes: dict[int, Node] = {}
    names: dict[int, str] = {}  # {minted id: user label} — for external trees; empty for ZOMBI ones
    written: dict[int, str] = {}  # {id: the label as written} — for ZOMBI trees; empty for external
    counter = 0

    def _mint(p: _P) -> int:
        nonlocal counter
        if all_labelled:
            m = _ZOMBI_LABEL.match(p.name)      # `all_labelled` is exactly "every label matched"
            assert m is not None
            return int(m.group(2))
        i_ = counter
        counter += 1
        return i_

    # first pass: assign ids (parents before children) and absolute times from durations
    def _build(p: _P, parent: int | None, birth: float) -> int:
        nid = _mint(p)
        if nid in nodes:
            raise ValueError(f"duplicate node id n{nid} in the Newick (labels must be unique)")
        end = birth + p.length
        kids = [_build(c, nid, end) for c in p.children]     # arity was checked while parsing
        child_ids = (kids[0], kids[1]) if kids else ()
        nodes[nid] = Node(nid, parent, birth, end, child_ids)  # fate filled in below
        if all_labelled:
            written[nid] = p.name              # the label as the file spells it: n<id> or e<id>
        elif p.name:                           # external tree: keep the user's label for the name-map
            names[nid] = p.name
        return nid

    root_id = _build(root_p, None, 0.0)

    # second pass: fates. Internal nodes are always speciations; the tips depend on the tree kind.
    for n in nodes.values():
        if n.children:
            n.fate = "speciation"
    leaves = [n for n in nodes.values() if not n.children]

    if assume_extant:
        # geometric callers (the tree transforms, treedist) don't use fate: take every tip as extant
        # and keep the branch lengths exactly as read — no ultrametric check, no snap. This is how a
        # non-ultrametric input (an inferred phylogram, a rounding-noisy dated tree) loads for them.
        for n in leaves:
            n.fate = "extant"
        return Tree(nodes, root_id), names

    # A tree with no time in it cannot carry a timed process. Every rate here is per unit time, so
    # on a topology-only tree each one fires zero times and the run "succeeds" having simulated
    # nothing — a genome per node, a gene tree per family, and an event log holding only the
    # originations and the speciations the topology itself forced. That is a hard failure to notice
    # from the outside, so it is refused at the door. The geometric callers above are exempt: a
    # topology with no lengths is a perfectly good input to a tree comparison.
    if sum(n.end_time - n.birth_time for n in nodes.values()) <= 0.0:
        raise ValueError(
            "this tree has no branch lengths, so it spans no time and nothing can evolve along it: "
            "every rate is per unit time, and on a tree of total length zero every one of them "
            "fires zero times. Give the tree branch lengths in time units — or simulate one with "
            "'zombi2 species'.")

    if all_labelled:
        if tip_fates is not None:
            # the run's own species_fates.tsv (or a --tip-fates file) states each tip's fate directly,
            # keyed by the same label the tree carries — both come from node_label, so they agree.
            # It is authoritative: it is the only thing that can tell an UNSAMPLED survivor from an
            # extant one (both sit at the present, and both are written n<id>).
            _assign_fates_from_map(leaves, {n.id: written[n.id] for n in leaves}, tip_fates,
                                   source="tip fates")
            return Tree(nodes, root_id), names
        # An e<id> label IS the fate: to_newick writes it for a lineage that died, so a ZOMBI tree
        # says so itself and needs neither the sibling fate table nor a guess from tip depth — which
        # is what makes the file survive being moved, copied or emailed. A tree written before the
        # letter existed carries none of them and falls through to the depth rule below, where it
        # behaves exactly as it did.
        if any(written[n.id][:1] == "e" for n in leaves):
            for n in leaves:
                n.fate = "extinct" if written[n.id][:1] == "e" else "extant"
            return Tree(nodes, root_id), names
        # no fate table: a tip is extinct if it ends before the present (the greatest end_time). The
        # tolerance is depth-relative — ``to_newick``'s default writes lengths exactly, but a tree
        # written at a fixed ``precision=`` (or hand-edited, or from a third party) accumulates
        # rounding along a root-to-tip path, so a tip at the present can fall a little short of the
        # max, far below any real extinction gap. The margin here is deliberately wide: it costs
        # nothing and absorbs all three. This cannot recover an unsampled tip (it sits at the
        # present, so it reads back extant) — pass the fate table for that.
        present = max(n.end_time for n in nodes.values())
        tol = max(1e-9, 1e-4 * present)
        for n in leaves:
            n.fate = "extinct" if n.end_time < present - tol else "extant"
        return Tree(nodes, root_id), names

    # an external tree: ultrametric ⟺ every tip is contemporaneous (all extant); otherwise the
    # differing depths could be extinctions or early samples, which we refuse to guess (SPEC decision).
    depths = [n.end_time for n in leaves]  # root sits at 0, so a tip's depth is its end_time
    height = max(depths)
    gap = max(depths) - min(depths)
    if gap <= max(1e-12, 1e-6 * height):  # ultrametric
        for n in leaves:
            n.fate = "extant"
        return Tree(nodes, root_id), names

    _assign_external_fates(leaves, names, tip_fates, gap)
    return Tree(nodes, root_id), names

zombi2.tree.with_stem

with_stem(tree: Tree, length: float, *, mode: str = 'set') -> Tree

Return a copy whose stem — the branch above the crown (root) — is length (mode="set") or is extended by length (mode="add"). Every other branch length is unchanged, so to_newick writes the new stem as )n<root>:<stem>; and nothing below moves.

Source code in zombi2/tree.py

def with_stem(tree: Tree, length: float, *, mode: str = "set") -> Tree:
    """Return a copy whose **stem** — the branch above the crown (root) — is ``length`` (``mode="set"``)
    or is extended by ``length`` (``mode="add"``). Every other branch length is unchanged, so
    ``to_newick`` writes the new stem as ``)n<root>:<stem>;`` and nothing below moves."""
    if not math.isfinite(length):
        raise ValueError(f"stem length must be finite, got {length!r}")
    if mode not in ("set", "add"):
        raise ValueError(f"mode must be 'set' or 'add', got {mode!r}")
    out = _copy(tree)
    root = out.nodes[out.root]
    root.birth_time = (root.end_time - length) if mode == "set" else (root.birth_time - length)
    if root.end_time - root.birth_time < 0:
        raise ValueError("resulting stem is negative")
    return out

zombi2.tree.make_ultrametric

make_ultrametric(tree: Tree, *, tol: float = 0.001) -> Tree

Return a copy in which every tip sits at the present (exactly ultrametric), by extending the terminal branches to a common depth. Snaps only when the tip-depth spread is within tol of the tree height — i.e. rounding; a larger spread raises, because differing tip depths then carry real signal (extinct lineages or serial samples) that this must not silently flatten.

A fixed precision= undoes this. The snap is exact here — tip depths agree to about 1e-16 — and to_newick's default (precision=None) writes every length back exactly, so the round trip keeps it. A fixed precision= does not: a depth is a sum of branch lengths, so to_newick(precision=7) reintroduces a spread of roughly 1e-6 on an ordinary tree, more than enough for ape::is.ultrametric() to reject the file this function was called to produce. zombi2 tools tree --round writes at the default, so the file it produces is still ultrametric.

Source code in zombi2/tree.py

def make_ultrametric(tree: Tree, *, tol: float = 1e-3) -> Tree:
    """Return a copy in which every tip sits at the present (exactly ultrametric), by extending the
    terminal branches to a common depth. Snaps only when the tip-depth spread is within ``tol`` of
    the tree height — i.e. rounding; a larger spread raises, because differing tip depths then carry
    real signal (extinct lineages or serial samples) that this must not silently flatten.

    **A fixed ``precision=`` undoes this.** The snap is exact here — tip depths agree
    to about 1e-16 — and ``to_newick``'s default (``precision=None``) writes every length back
    exactly, so the round trip keeps it. A *fixed* ``precision=`` does not: a depth is a *sum* of
    branch lengths, so ``to_newick(precision=7)`` reintroduces a spread of roughly 1e-6 on an
    ordinary tree, more than enough for ``ape::is.ultrametric()`` to reject the file this function was
    called to produce. ``zombi2 tools tree --round`` writes at the default, so the file it produces
    is still ultrametric."""
    depth = _depths(tree)
    tips = [i for i, n in tree.nodes.items() if not n.children]
    lo, hi = min(depth[i] for i in tips), max(depth[i] for i in tips)
    if hi > 0 and (hi - lo) > tol * hi:
        raise ValueError(
            f"tip depths differ by {hi - lo:.3g} (> {tol:g} × height {hi:.3g}); this is more than "
            "rounding — the tips are not contemporaneous (extinct lineages or serial samples), so "
            "there is no ultrametric tree to snap to")
    out = _copy(tree)
    for i in tips:
        nd = out.nodes[i]
        parent_depth = depth[i] - (nd.end_time - nd.birth_time)   # = depth of this tip's parent
        nd.end_time = nd.birth_time + (hi - parent_depth)         # so the tip lands at depth hi
    return out

zombi2.tree.rescale

rescale(tree: Tree, *, height: float | None = None, factor: float | None = None) -> Tree

Return a copy with every branch length scaled — either so the root-to-tip height equals height, or by a raw factor. Exactly one of the two must be given.

Source code in zombi2/tree.py

def rescale(tree: Tree, *, height: float | None = None, factor: float | None = None) -> Tree:
    """Return a copy with every branch length scaled — either so the root-to-tip height equals
    ``height``, or by a raw ``factor``. Exactly one of the two must be given."""
    if (height is None) == (factor is None):
        raise ValueError("pass exactly one of height= or factor=")
    if factor is None:
        assert height is not None                # exactly one of the two was given, checked above
        depth = _depths(tree)
        current = max(depth[i] for i, n in tree.nodes.items() if not n.children)
        if current <= 0:
            raise ValueError("tree has zero height; cannot scale it to a target height")
        factor = height / current
    if factor < 0:
        raise ValueError(f"scale factor must be non-negative, got {factor}")
    out = _copy(tree)
    for nd in out.nodes.values():
        nd.birth_time *= factor
        nd.end_time *= factor
    return out

zombi2.tree.relative_evolutionary_divergence

relative_evolutionary_divergence(tree: Tree) -> dict[int, float]

Relative Evolutionary Divergence (Parks et al. 2018) of every node — root 0.0, leaves 1.0, keyed by node id. Walking root-outward, a node sits at RED(parent) + a/(a+b)·(1 − RED(parent)) where a is its branch and b the mean branch-length distance from it to the leaves of its subtree. RED is invariant to a global rescaling, so a rate-distorted phylogram reads as an approximate relative timeline; on an ultrametric tree it returns each node's exact relative age. A zero-length branch passes the parent's value straight down.

Source code in zombi2/tree.py

def relative_evolutionary_divergence(tree: Tree) -> dict[int, float]:
    """Relative Evolutionary Divergence (Parks et al. 2018) of every node — root ``0.0``, leaves
    ``1.0``, keyed by node id. Walking root-outward, a node sits at ``RED(parent) + a/(a+b)·(1 −
    RED(parent))`` where ``a`` is its branch and ``b`` the mean branch-length distance from it to the
    leaves of its subtree. RED is invariant to a global rescaling, so a rate-distorted phylogram reads
    as an approximate relative timeline; on an ultrametric tree it returns each node's exact relative
    age. A zero-length branch passes the parent's value straight down."""
    nodes = tree.nodes
    order = _preorder(tree)
    if not order:
        raise ValueError("empty tree — nothing to compute RED on")

    def length(i: int) -> float:
        a = nodes[i].end_time - nodes[i].birth_time
        if a < 0.0:
            raise ValueError(f"negative branch length ({a}) above node n{i}")
        return a

    mean_tip_dist: dict[int, float] = {}
    n_leaves: dict[int, int] = {}
    for i in reversed(order):                       # child before parent
        kids = nodes[i].children
        if not kids:
            mean_tip_dist[i] = 0.0
            n_leaves[i] = 1
            continue
        total = 0.0
        k = 0
        for c in kids:
            total += n_leaves[c] * (length(c) + mean_tip_dist[c])
            k += n_leaves[c]
        mean_tip_dist[i] = total / k
        n_leaves[i] = k

    red: dict[int, float] = {}
    for i in order:                                 # parent before child
        p = nodes[i].parent
        if p is None:
            red[i] = 0.0
            continue
        a, b, pr = length(i), mean_tip_dist[i], red[p]
        red[i] = pr + (a / (a + b)) * (1.0 - pr) if (a + b) > 0.0 else pr
    return red

zombi2.tree.red_scaled

red_scaled(tree: Tree) -> Tree

Return a copy whose node depths are their RED — ultrametric on [0, 1], root at 0, every tip at 1. Branch lengths become RED increments. This is the tree GTDB-style rank normalisation reads (relative_evolutionary_divergence() gives the raw per-node values).

Source code in zombi2/tree.py

def red_scaled(tree: Tree) -> Tree:
    """Return a copy whose node depths **are** their RED — ultrametric on ``[0, 1]``, root at 0, every
    tip at 1. Branch lengths become RED increments. This is the tree GTDB-style rank normalisation
    reads (`relative_evolutionary_divergence()` gives the raw per-node values)."""
    red = relative_evolutionary_divergence(tree)
    out = _copy(tree)
    for i, nd in out.nodes.items():
        nd.birth_time = 0.0 if nd.parent is None else red[nd.parent]
        nd.end_time = red[i]
    return out

zombi2.tree.gamma_statistic

gamma_statistic(tree: Tree) -> float

Pybus & Harvey's γ — where a dated tree's branching times sit relative to what a constant rate would give.

Standard normal under constant-rate pure birth, so a value near 0 is the null. It goes negative when speciation slows toward the present, because the branching times then bunch up early, and positive when it accelerates. The tree must be dated and ultrametric: γ reads the waiting times between splits, so branch lengths in substitutions mean nothing here. Extinct lineages must be pruned first — γ is defined on the reconstructed tree.

Needs at least four tips: the statistic divides by n - 2 and by the tree's total branch length, neither of which is usable below that.

Source code in zombi2/tree.py

def gamma_statistic(tree: Tree) -> float:
    """Pybus & Harvey's γ — where a dated tree's branching times sit relative to what a constant rate
    would give.

    Standard normal under constant-rate pure birth, so a value near 0 is the null. It goes negative
    when speciation slows toward the present, because the branching times then bunch up early, and
    positive when it accelerates. The tree must be **dated and ultrametric**: γ reads the waiting
    times between splits, so branch lengths in substitutions mean nothing here. Extinct lineages must
    be pruned first — γ is defined on the reconstructed tree.

    Needs at least four tips: the statistic divides by ``n - 2`` and by the tree's total branch
    length, neither of which is usable below that.
    """
    splits = sorted(nd.end_time for nd in tree.nodes.values() if nd.children)
    leaves = [nd.end_time for nd in tree.nodes.values() if not nd.children]
    if not leaves:
        raise ValueError("empty tree — nothing to compute gamma on")
    n = len(splits) + 1
    if n < 4:
        raise ValueError(f"gamma needs at least 4 tips, this tree has {n}")
    present = max(leaves)
    spread = present - min(leaves)
    if spread > 1e-6 * max(present, 1.0):
        raise ValueError("gamma is defined on a dated ultrametric tree, and this one's tips differ "
                         f"in depth by {spread:g} — prune the extinct lineages, or pass a dated tree "
                         "rather than a phylogram")
    # inter[k] is the waiting time during which k + 2 lineages were alive
    inter = [splits[j] - splits[j - 1] for j in range(1, n - 1)] + [present - splits[-1]]
    partial, total = [], 0.0
    for k, g in enumerate(inter):
        total += (k + 2) * g
        partial.append(total)
    if total <= 0.0:
        raise ValueError("gamma needs a tree with positive branch lengths")
    return (sum(partial[:-1]) / (n - 2) - total / 2) / (total * (1 / (12 * (n - 2))) ** 0.5)

zombi2.tree.distance

distance(a: Tree, b: Tree, *, metric: str = 'rf') -> float

Distance between two rooted trees over their shared tips (matched by node id). Raises if the two leaf sets differ. metric: "rf" (Robinson–Foulds — the number of clades in one tree but not the other), "rf-normalized" (that count over the total number of non-trivial clades), or "branch-score" (Kuhner–Felsenstein — √Σ(branch-length difference)² over all clades, terminal branches included).

Source code in zombi2/tree.py

def distance(a: Tree, b: Tree, *, metric: str = "rf") -> float:
    """Distance between two **rooted** trees over their shared tips (matched by node id). Raises if the
    two leaf sets differ. ``metric``: ``"rf"`` (Robinson–Foulds — the number of clades in one tree but
    not the other), ``"rf-normalized"`` (that count over the total number of non-trivial clades), or
    ``"branch-score"`` (Kuhner–Felsenstein — √Σ(branch-length difference)² over all clades, terminal
    branches included)."""
    la = frozenset(i for i, n in a.nodes.items() if not n.children)
    lb = frozenset(i for i, n in b.nodes.items() if not n.children)
    if la != lb:
        raise ValueError(
            f"the two trees have different leaf sets ({len(la)} vs {len(lb)} tips, {len(la ^ lb)} not "
            "shared) — treedist needs the same taxa, identically labelled, on both trees")
    ca, cb = _clades(a), _clades(b)
    n = len(la)
    if metric in ("rf", "rf-normalized"):
        sa = {c for c in ca if 2 <= len(c) <= n - 1}
        sb = {c for c in cb if 2 <= len(c) <= n - 1}
        rf = len(sa ^ sb)
        if metric == "rf":
            return float(rf)
        denom = len(sa) + len(sb)
        return float(rf / denom) if denom else 0.0
    if metric == "branch-score":
        return float(sum((ca.get(c, 0.0) - cb.get(c, 0.0)) ** 2 for c in set(ca) | set(cb)) ** 0.5)
    raise ValueError(f"unknown metric {metric!r}; choose 'rf', 'rf-normalized', or 'branch-score'")