zombi2.genomes

Level 2: genomes evolving along the species tree. The level has three resolutions — family ⊂ ordered ⊂ nucleotide — and one entry point each; more detail costs more compute, so the resolution is a dial.

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, families=None, max_family_size=10, joint=False, seed=None, parallel=False, stream_to=None, outputs=None, compact=False, bundle=None, progress=False, **retired) -> '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 Tree, or a 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 Recipients().weighted_by(driver, mapping) (weighted by an evolved value; 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 origin. families=[family("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 scaled_by("genomes:toxin", …) reads. Deterministic given seed.

A family with rates of its own. families=[family("IS1", transfer=PerCopy(1.5), loss=0.02)] declares a named family and gives it its own rates, so one family can behave nothing like the rest of the genome — a mobile element that transfers constantly, a core gene that is almost never lost. Whatever a family leaves out falls back to the run's rate for that event, so families=[family("toxin")] is a family declared for its name alone and is exactly families=[family("toxin")]. family() also takes origin= and module=, which are origins= and modules= said on the family itself. Origination takes no per-family value: when it is read the family does not exist yet to have one.

A written rate is that family's rate, in the rate's own units. The run's rate does not reach it, and neither does a varying_among draw meant to vary the run's rate — not the per-family one, and not the per-lineage one either, so a family that writes its own loss keeps that loss on a branch every other family loses fast on (SPEC §5's argument for set_by, one level down). Two things are refused for now: a family's rate carrying a verb of its own, and a family's rate beside a driven run rate; both arrive with the joint step. So does the per-family engine — parallel= and stream_to= build one set of rates for the whole run, so they refuse a family that writes its own.

A family placed by hand. origins=[("n5", 0.4)] originates a family on lineage n5 at time 0.4 — the same event the origination rate produces, at a point you choose rather than one that is drawn. It adds to the run: whatever initial_families and origination give you is still there, so origins with both of those at 0 is a tree carrying exactly the families you placed, and origins beside them is an ordinary genome with one family started where you want it. The time is the run's own clock and must fall inside that lineage's life (None puts it at the branch's start); a placed family is an ordinary family from that instant on, so it duplicates, transfers and is lost like any other, and gets its gene tree the same way. The ids come straight after the initial and named ones, in the order you wrote the origins, on either engine — so the family you placed is initial_families + len(family_names). See resolve_origins.

Conditioning (a trait drives a rate). Any of the four rates may be driven by another level — loss = PerCopy(0.25).scaled_by("trait_events.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=("events",)), which writes trait_events.tsv). 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 = Recipients().weighted_by(driver, 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, a weight, not a rate). 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 _do_transfer()). The two driven arguments are independent and may be used together or apart.

In a joint=True run the driver can also be gene content the run is building: Recipients().weighted_by("genomes:A", {"present": 20.0, "absent": 1.0}) weighs each candidate by whether it carries A when the transfer fires. A declared family can give its own transfer_to (see family()), which is used for its copies instead of the run's.

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 a weighted transfer_to runs on the per-family engine too — conditioning does not couple families, so nothing here forces a fallback. 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 StreamedRun handle comes back instead of a FamilyGenomesResult. outputs= picks which files, as FamilyGenomesResult.write() takes them minus summary (a streamed run writes no genome_summary.json); the default is all six. It is the per-family engine, and outputs without stream_to is an error. compact=True writes the gene trees as one gene_trees.tsv, as write(compact=True) does; without stream_to it is an error too.

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, families=None, max_family_size=10, joint=False,
                            seed=None, parallel=False, stream_to=None, outputs=None,
                            compact=False, bundle=None, progress=False,
                            **retired) -> "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 `Tree`, or a
    `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 ``Recipients().weighted_by(driver, mapping)`` (weighted by an evolved value; 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 origin. ``families=[family("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
    ``scaled_by("genomes:toxin", …)`` reads. Deterministic given ``seed``.

    **A family with rates of its own.** ``families=[family("IS1", transfer=PerCopy(1.5), loss=0.02)]``
    declares a named family and gives it its own rates, so one family can behave nothing like the rest
    of the genome — a mobile element that transfers constantly, a core gene that is almost never lost.
    Whatever a family leaves out falls back to the run's rate for that event, so
    ``families=[family("toxin")]`` is a family declared for its name alone and is exactly
    ``families=[family("toxin")]``. `family()` also takes ``origin=`` and ``module=``, which are
    ``origins=`` and ``modules=`` said on the family itself. Origination takes no per-family value:
    when it is read the family does not exist yet to have one.

    A written rate **is** that family's rate, in the rate's own units. The run's rate does not reach
    it, and neither does a ``varying_among`` draw meant to vary the run's rate — not the per-family
    one, and not the per-lineage one either, so a family that writes its own loss keeps that loss on
    a branch every other family loses fast on (SPEC §5's argument for ``set_by``, one level down). Two things are refused for now: a
    family's rate carrying a verb of its own, and a family's rate beside a *driven* run rate; both
    arrive with the joint step. So does the per-family engine — ``parallel=`` and ``stream_to=``
    build one set of rates for the whole run, so they refuse a family that writes its own.

    **A family placed by hand.** ``origins=[("n5", 0.4)]`` originates a family on lineage ``n5`` at
    time ``0.4`` — the same event the ``origination`` rate produces, at a point you choose rather
    than one that is drawn. It **adds to** the run: whatever ``initial_families`` and ``origination``
    give you is still there, so ``origins`` with both of those at 0 is a tree carrying exactly the
    families you placed, and ``origins`` beside them is an ordinary genome with one family started
    where you want it. The time is the run's own clock and must fall inside that lineage's life
    (``None`` puts it at the branch's start); a placed family is an ordinary family from that instant
    on, so it duplicates, transfers and is lost like any other, and gets its gene tree the same way.
    The ids come straight after the initial and named ones, in the order you wrote the origins, on
    either engine — so the family you placed is ``initial_families + len(family_names)``. See
    `resolve_origins`.

    **Conditioning (a trait drives a rate).** Any of the four rates may be *driven by another level* —
    ``loss = PerCopy(0.25).scaled_by("trait_events.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=("events",))``, which writes
    ``trait_events.tsv``). 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 =
    Recipients().weighted_by(driver, 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, a weight, not a rate). 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 `_do_transfer()`). The two driven arguments are
    independent and may be used together or apart.

    In a ``joint=True`` run the driver can also be gene content the run is building:
    ``Recipients().weighted_by("genomes:A", {"present": 20.0, "absent": 1.0})`` weighs each candidate
    by whether it carries ``A`` when the transfer fires. A declared family can give its own
    ``transfer_to`` (see `family()`), which is used for its copies instead of the run's.

    ``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 a weighted ``transfer_to`` runs on the per-family engine too — conditioning does not couple
    families, so nothing here forces a fallback. 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 `StreamedRun` handle comes
    back instead of a ``FamilyGenomesResult``. ``outputs=`` picks which files, as
    `FamilyGenomesResult.write()` takes them minus ``summary`` (a streamed run writes no
    ``genome_summary.json``); the default is all six. It is the per-family engine, and
    ``outputs`` without ``stream_to`` is an error. ``compact=True`` writes the gene trees as one
    ``gene_trees.tsv``, as ``write(compact=True)`` does; without ``stream_to`` it is an error too.
    """
    # First line of the call, before any work: a worker re-importing an unguarded script would
    # otherwise repeat the whole run before dying at its own pool. See `refuse_worker_reentry`.
    refuse_worker_reentry(parallel)
    compact = compact_or_bundle(compact, bundle, "simulate_genomes_family")   # its 0.52.0 name
    tree = as_tree(tree, level="genomes")
    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)
    # A scope says *per what* a rate is counted, and the three copy-consuming events implement two
    # answers. **Per copy** (the default) puts every copy independently at risk, so a bigger genome
    # turns over faster. **Per lineage** is a fixed budget: the lineage duplicates or loses at its
    # rate whatever its genome size — and then the total and the pick both have to change, or the
    # rate would say one thing and the picking another. Origination is per lineage alone, because it
    # is the rate at which families are CREATED: per copy it would be base × 0 in an empty genome, a
    # silent no-op. A driver is read per lineage on all four events; on transfer the driven
    # lineage is the DONOR (who receives is the separate transfer_to choice, below).
    for label, rate, legal in (("duplication", dup, (PerCopy, PerLineage)),
                               ("transfer", tra, (PerCopy, PerLineage)),
                               ("loss", los, (PerCopy, PerLineage)),
                               ("origination", org, (PerLineage,))):
        # `rate.scope` holds the scope CLASS, not an instance — a scope constructor returns the rate
        # itself — so this is an identity test against the legal set rather than an isinstance one.
        assert rate.scope is not None            # as_rate fills the level's default where none was written
        if rate.scope not in legal:
            raise ValueError(
                f"{label} has a {rate.scope.__name__} scope, but the family genome engine "
                f"takes {' or '.join(s.__name__ for s in legal)} for {label}."
            )
        for m in rate.modifiers:
            if m.reads == (DRAWN, "families") and label == "origination":
                raise ValueError(
                    "origination carries a per-family draw, 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. "
                    "Write varying_among('families', …) on duplication, transfer or loss; writing one "
                    "such object on several of them gives a family-wide tempo, since one object is "
                    "one draw.")
            if isinstance(m, Driven):
                check_not_a_kernel(m.mapping, label=label)
            if is_implemented(m, IMPLEMENTED_MODIFIERS, "genomes.family"):
                continue
            raise ValueError(
                f"{label} carries {describe(m)}, which the family genome engine does not "
                f"support. It takes changing_at (skyline), scaled_by (a conditioned or joint driver), "
                f"set_by (a driver that replaces the base), varying_among('families', …) "
                f"(per-family heterogeneity) and varying_among('lineages', …) (one multiplier per "
                f"species branch, drawn or inherited)."
            )
    for label, rate in (("duplication", dup), ("transfer", tra), ("loss", los),
                        ("origination", org)):
        rate.check_one_base(label)
        # 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 branch's factor to be independent of
        # its parent's and inherited from it at once — there is no model there to implement.
        check_one_memory(tuple(m for m, _ in rate.carried_modifiers(unit="lineages")),
                         label=label, unit="lineages")
    # A per-family draw under a per-lineage scope is refused because what it should MEAN is a
    # modelling decision nobody has taken, not because it is hard. Under PerCopy a family's
    # multiplier scales each copy's own rate, so a genome full of fast families turns over faster —
    # the multiplier moves the total. Under PerLineage the total is fixed by definition, so the same
    # multiplier could only decide *which* copy the event takes, normalised within the lineage. Those
    # are two different models, and running one while the user wrote the other is exactly the silent
    # mismatch this engine refuses everywhere else.
    # The check is over the whole RUN, not per rate, and getting that wrong was a real bug: a
    # per-family draw anywhere makes the engine take its per-family path for **every** gene rate,
    # summing each one over the live copies — so a `PerLineage` rate elsewhere in the same run was
    # silently counted per copy while nothing on the page said so. One draw on duplication was
    # enough to turn a `PerLineage` loss back into a `PerCopy` one.
    per_lineage_here = [lbl for lbl, r in (("duplication", dup), ("transfer", tra), ("loss", los))
                        if r.scope is PerLineage]
    drawn_here = [lbl for lbl, r in (("duplication", dup), ("transfer", tra), ("loss", los))
                  if any(m.reads == (DRAWN, "families") for m in r.modifiers)]
    if per_lineage_here and drawn_here:
        raise ValueError(
            f"{', '.join(per_lineage_here)} is PerLineage while {', '.join(drawn_here)} carries a "
            f"per-family draw, and the two cannot share a run. Under PerCopy a family's multiplier "
            f"scales each copy's rate, so it changes the lineage's total; under PerLineage the total "
            f"is fixed whatever the genome holds, so the multiplier could only choose which copy is "
            f"taken. Those are different models and the choice is not made yet — write PerCopy "
            f"throughout for the first, or drop the per-family draw for the second.")
    # the choice (SPEC §5), validated in the one place all three resolutions share: the mapping's
    # numbers are weights over the candidate recipients, never a rate multiplier
    transfer_to = resolve_transfer_to(transfer_to)
    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}")
    check_no_retired_keywords(retired, where="simulate_genomes_family")
    # every named family, from the one list that declares them
    declared, module_map, started_named = resolve_families(families, tree)
    family_names = [f.name for f in declared]
    # a family's own recipient rule, checked by the same resolver as the run's transfer_to
    fam_transfer_to = resolve_family_transfer_to(declared)
    # what each family sets for itself; empty unless some family writes a rate, and then the engine
    # takes the path that sums those beside the run's (see `_family_weights`)
    # (a rate carrying a verb is kept apart, in `fam_driven`, and read per lineage — see `own_sums`)
    fam_own, fam_driven = resolve_family_rates(
        declared, {"duplication": dup, "transfer": tra, "loss": los})
    any_written = bool(fam_own) or bool(fam_driven)
    if any_written:
        for key, rate in (("duplication", dup), ("transfer", tra), ("loss", los)):
            if key not in fam_own and key not in fam_driven:
                continue
            assert rate.scope is not None        # `as_rate` filled the level's default above
            if rate.scope is not PerCopy:
                raise ValueError(
                    f"a family writes its own {key}, but the run's {key} is {rate.scope.__name__}. "
                    f"The two are summed over the same copies, so both are counted per copy — write "
                    f"PerCopy for the run's {key}, or drop the family's.")
            # a driven run rate composes: the driver scales the run's rate, and a family that states
            # its own has no run rate in it to scale (see `_driven_weights`)

    # 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)

    # conditioning: a rate written with scaled_by 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 (a str driver) or an object handed over in memory (a trait result, or a genome's presence /
    # completion). No driven rate ⇒ this is empty
    # and the loop stays byte-identical to an undriven run.
    dup_mods, los_mods = _driven_mods(dup), _driven_mods(los)
    org_mods, tra_mods = _driven_mods(org), _driven_mods(tra)
    # A family's own rate reads its drivers through the same checks as the run's rates: a family it
    # names must be declared, a live name needs joint=True, and a finished driver becomes a trajectory
    # the loop steps at. Appended last, so a run with no such rate resolves its drivers as before.
    fam_rate_mods = [m for table in fam_driven.values() for rate in table.values()
                     for m in _driven_mods(rate)]
    all_mods = (*dup_mods, *los_mods, *org_mods, *tra_mods, *fam_rate_mods)
    # Two kinds of driver, told apart by what the driver *is* (SPEC §5). A finished one — a file, a
    # grown result — was produced before this run started, and the run is conditioned. A **live**
    # name is gene content this run is itself producing, and the run is joint: there is no order to
    # simulate the two in, because they are the same run. Only the first can be resolved to a
    # trajectory; the second is read off the live genome as the loop goes.
    live_mods = [m for m in all_mods if names_a_live_level(m.driver)]
    file_mods = [m for m in all_mods if not names_a_live_level(m.driver)]
    # A transfer_to weighted by live gene content is read when a transfer fires, off the genomes the
    # run is building. It is checked with the rates' live drivers and makes the run joint too, but it
    # stays out of `live_keys`, because it moves no rate (see `prepare_transfer_to`).
    live_choices = [r for r in (transfer_to, *fam_transfer_to.values())
                    if isinstance(r, Driven) and names_a_live_level(r.driver)]
    live_keys = resolve_live_drivers(
        live_mods, set(family_names), joint=joint, choice_mods=live_choices,
        modules={name: len(members) for name, members in (module_map or {}).items()})
    # driver key → its Driven (deduped, so a driver shared across rates resolves once);
    # the modifier rather than the driver itself, because the driver's step rides on the modifier
    by_key: dict[object, "Driven"] = {}
    for m in file_mods:
        by_key.setdefault(m.key, m)
    resolved = {}
    if by_key:
        resolved = {key: resolve_driver(m.driver, tree, step=m.step, path=m.path, level="genomes.family")
                    for key, m 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 undriven model — refuse it here,
        # naming the driver, rather than let it pass as a driven run
        for m in file_mods:
            label = m.driver if isinstance(m.driver, str) else f"<{type(m.driver).__name__}>"
            check_mapping_fires(m.mapping, resolved[m.key].states(), driver_label=label)
    # `trajs` is the drivers that move a RATE: they alone make the loop per-lineage and set the
    # Gillespie horizon. It is built BEFORE transfer_to is prepared, and that order is
    # load-bearing — a driven transfer_to changes no rate, so its trajectory must not end up here
    # adding horizon breakpoints (see prepare_transfer_to). `resolved` is passed along as the driver
    # cache, so a driver shared between a rate and transfer_to is loaded once.
    trajs = dict(resolved)

    def prepare_choice(rule):
        """``(groups, trajectory)`` for a recipient rule. A rule reading live gene content has nothing
        to prepare here: its reader is built once the family ids and the genomes exist (below)."""
        if isinstance(rule, Driven) and names_a_live_level(rule.driver):
            return None, None
        return prepare_transfer_to(tree, rule, resolved, level="genomes.family")

    group_of, to_traj = prepare_choice(transfer_to)
    fam_choice_prepared = {i: (rule, *prepare_choice(rule)) for i, rule in fam_transfer_to.items()}

    # 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). The
    # drivers and clade painting above are resolved once, before the split, and handed to whichever
    # engine runs — they are shared validation and shared input, not one engine's business. A
    # configuration neither engine covers still returns None there, so this serial reference loop runs
    # unchanged (decision A); a streamed run raises instead, being unable to fall back without pulling
    # the whole thing into memory.
    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 compact and stream_to is None:
        raise ValueError(
            "compact applies to a streamed run (stream_to=DIR), which writes the files itself; for "
            "an in-memory run pass it when you call result.write(compact=True).")
    seed = resolve_seed(seed)     # drawn if none was given, so either engine below records it
    if (parallel or stream_to is not None) and (live_keys or live_choices):
        # The one thing that engine's whole design rests on: a family's history depends on no other
        # family, so each can be evolved alone. A rate reading the genome's own content is exactly
        # that dependence, so this is a refusal rather than a fallback.
        raise ValueError(
            "a joint genome run cannot use the per-family engine (parallel= / stream_to=), which "
            "evolves each family in its own process because families do not affect each other. A "
            "rate or transfer_to reading live gene content is that effect. Drop parallel / stream_to.")
    if (parallel or stream_to is not None) and started_named:
        # Pass 1 of that engine enumerates every family's origination up front, seeding the declared
        # ones at the root; a family that arrives partway down is not in that enumeration, and
        # threading it through would renumber the families the serial engine mints. Stated rather
        # than silently dropped.
        raise ValueError(
            "family(origin=...) does not run on the per-family engine (parallel= / stream_to=), "
            "which enumerates every family's origination before it starts. Drop parallel / "
            "stream_to, or let the origination rate place the family.")
    if (parallel or stream_to is not None) and any_written:
        # not a fallback: that engine evolves one family per process against a context built once for
        # the whole run, so a per-family rate has nowhere to live in it. Saying so beats silently
        # running every family at the run's rate.
        raise ValueError(
            "a family writing its own rate cannot run on the per-family engine (parallel= / "
            "stream_to=), which builds one set of rates for the whole run and evolves each family "
            "against it. Drop parallel / stream_to, or give every family the run's rate.")
    if (parallel or stream_to is not None) and fam_transfer_to:
        raise ValueError(
            "a family with its own transfer_to cannot run on the per-family engine (parallel= / "
            "stream_to=), which builds one recipient rule for the whole run and evolves each family "
            "against it. Drop parallel / stream_to, or give every family the run's transfer_to.")
    # What each rate carries among lineages: a factor per species branch, shared by every family
    # passing through it — the other axis from the per-family draws (`zombi2.genomes.multipliers`).
    # Read here, above the split, because both engines need it; each draws its own numbers.
    lin_by = {"duplication": tuple(m for m, _ in dup.carried_modifiers(unit="lineages")),
              "transfer": tuple(m for m, _ in tra.carried_modifiers(unit="lineages")),
              "loss": tuple(m for m, _ in los.carried_modifiers(unit="lineages")),
              "origination": tuple(m for m, _ in org.carried_modifiers(unit="lineages"))}
    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, placed=[], modules=module_map, cap=cap,
            seed=seed, parallel=parallel,
            progress=progress, stream_to=stream_to, outputs=outputs, compact=compact,
            trajs=trajs, to_traj=to_traj, group_of=group_of,
            driven={"duplication": bool(dup_mods), "transfer": bool(tra_mods),
                    "loss": bool(los_mods), "origination": bool(org_mods)},
            lin_by=lin_by)
        if result is not None:
            return result

    rng, seed = stream("genomes", seed)     # the genomes level's own stream

    # The per-lineage multipliers this engine draws: one per species branch, drawn before any family
    # is minted, because they come from the tree, which exists before the run does. A run carrying
    # none draws nothing, so it is bit-identical to one from before this existed. On transfer the
    # branch is the DONOR's, as a driven transfer's is.
    lin_mult = draw_lineage_multipliers(lin_by, tree, rng)
    any_lineage = bool(lin_mult)

    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 created and then fixed for its whole life.
    # Whether a family's rates move together is decided by what was written: one `Random` object read
    # by two rates is one draw for both, two objects are two draws. Empty unless some rate carries
    # one, and then the engine takes its weighted path; a run carrying none draws nothing here.
    fam_by = {"duplication": tuple(m for m, _ in dup.carried_modifiers(unit="families")),
              "transfer": tuple(m for m, _ in tra.carried_modifiers(unit="families")),
              "loss": tuple(m for m, _ in los.carried_modifiers(unit="families"))}
    any_family = any(fam_by.values()) or any_written
    # A rate carrying nothing per family holds 1.0 for every family, so all such rates share one
    # empty table rather than each filling its own — which is what lets _FamilyWeights sum them once.
    # Sharing is off once some family writes a rate: a family may write its loss and not its
    # duplication, so the two tables then hold different numbers for the same family.
    no_variation: dict[int, float] = {}
    fam_mult: dict[str, dict[int, float]] = {
        key: ({} if (mods or any_written) else no_variation) for key, mods in fam_by.items()}
    #: per event, the families that set a rate of their own — their own per-copy rate, and 0.0 for
    #: every other family, which contributes through `fam_mult` instead. `None` when nobody does,
    #: and then every expression below is the one it always was (see `_family_weights`).
    fam_fixed: "dict[str, dict[int, float]] | None" = (
        {key: {} for key in fam_by} if any_written else None)
    #: per event, the families whose own rate carries a verb: minted id → that `Rate`. Such a family
    #: holds 0.0 in `fam_fixed` and in `fam_mult`, and its rate is read on a lineage when it is needed.
    fam_driven_ids: dict[str, dict[int, object]] = {key: {} for key in fam_by}
    #: per event, the families whose own rate, fixed or driven, replaces the run's: no draw reaches
    #: them there, so `family_multipliers` holds no number for them
    own_cells: dict[str, set[int]] = {key: set() for key in fam_by}
    any_family_driven = bool(fam_driven)

    def family_rate_on(rate, k: int) -> float:
        """A driven family rate read on lineage ``k`` at the current time: one copy's rate."""
        if any_driven:
            return rate.effective(copies=1, lineages=1, time=t, drivers=drivers[k])
        return rate.effective(copies=1, lineages=1, time=t)

    def own_rates(key: str, k: int):
        """The rate each family set for itself on this event, as lineage ``k`` reads it, or ``None``:
        what a copy pick reads, against the per-lineage *sums* that `own_sums` reads. Two shapes of the
        same information: a rate per family here, and per lineage the total over its live copies there.
        A driven family's rate differs between lineages, so for a driven family present on ``k``
        the value read on ``k`` replaces the 0.0 the table holds for it."""
        if fam_fixed is None:
            return None
        table = fam_driven_ids[key]
        if not table:
            return fam_fixed[key]
        return collections.ChainMap({f: family_rate_on(rate, k) for f, rate in table.items()
                                     if counts.holds(k, f)}, fam_fixed[key])

    def new_family(declared_at: "int | None" = None) -> int:
        """Mint a family id. ``declared_at`` is its index in ``declared`` for a named family, which is
        how it finds the rates it wrote; an anonymous family passes nothing and runs at the run's."""
        nonlocal family_counter
        f = family_counter
        family_counter += 1
        if any_family:
            # one draw per distinct modifier *object* for this family, shared across its rates: the
            # same `Random` object written on duplication and on loss means one number, so a fast
            # family is fast at both. Two separately built ones are two draws even with the same law.
            shared: dict[int, float] = {}
            for key, mods in fam_by.items():
                own = None if declared_at is None else fam_own.get(key, {}).get(declared_at)
                if own is not None:
                    # this family's rate IS the number it wrote, so the run's rate does not reach it
                    # and neither does a draw meant to vary the run's rate among families
                    assert fam_fixed is not None     # a written rate is exactly when it was built
                    fam_mult[key][f] = 0.0
                    fam_fixed[key][f] = own
                    own_cells[key].add(f)
                    continue
                driven_rate = None if declared_at is None else fam_driven.get(key, {}).get(declared_at)
                if driven_rate is not None:
                    # the same, except the number depends on the lineage and the time, so it is read
                    # where it is used (`own_sums`, `own_rates`) rather than stored here
                    assert fam_fixed is not None
                    fam_mult[key][f] = 0.0
                    fam_fixed[key][f] = 0.0
                    fam_driven_ids[key][f] = driven_rate
                    own_cells[key].add(f)
                    continue
                fam_mult[key][f] = math.prod(values_at_birth(mods, rng, shared))
                if fam_fixed is not None:
                    fam_fixed[key][f] = 0.0
        return f

    # Which of the three copy-consuming rates is a fixed per-lineage budget rather than a per-copy
    # risk. Read once: it decides both how the total is counted and how the affected lineage is
    # picked, and those two must never disagree.
    dup_per_lineage = dup.scope is PerLineage
    los_per_lineage = los.scope is PerLineage
    tra_per_lineage = tra.scope is PerLineage
    any_per_lineage = dup_per_lineage or los_per_lineage or tra_per_lineage

    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[GeneCopy]] = []
    pos: dict[int, int] = {}
    genomes: dict[int, tuple[GeneCopy, ...]] = {}
    events: list[GeneEdge] = []
    enter(alive, gen, pos, root.id, [])
    for _ in range(initial_families):  # lay down the origin's genome 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)
    named_starts: list[tuple[float, int, int]] = []
    for i, spec in enumerate(declared):
        fid = new_family(i)
        named[spec.name] = fid
        if i in started_named:
            # a declared family given an `origin` starts there rather than being seeded at the origin,
            # which is what `origins=` does for an anonymous one — the same event, with a name on it
            t_p, lineage = started_named[i]
            named_starts.append((t_p, lineage, fid))
            continue
        c = new_copy(fid)
        gen[0].append(c)
        events.append(GeneEdge(t, "origination", root.id, fid, c.id))
    # Such a family is not seeded into a genome: it arrives at its own time, in the loop below.
    starts = sorted(named_starts)
    start_i = 0                                                              # walked in time order
    total_copies = len(gen[0])
    initial_genome = tuple(gen[0])   # the run's starting genome: a snapshot before the stem runs

    any_driven = bool(trajs) or bool(live_keys)

    # Each live driver paired with what it reads, resolved once here: `live_keys` and `named` are both
    # fixed for the whole run, so the read below is a lookup rather than the same name taken apart
    # again on every one of the millions of reads a joint run makes.
    live_reads: list[tuple[str, int | tuple[int, ...] | None]] = [
        (src, live_target(src, named, module_map or {})) for src in live_keys]

    def live_value(fid: "int | tuple[int, ...] | None", k: int):
        """What a live driver reads on lineage ``k`` **right now** — the joint half of the driver
        mechanism (SPEC §2). A finished driver answers from a trajectory built before the run; this
        one answers from the genome the run is building. ``fid`` is the family the driver names, the
        ids of a module's families, or None for the driver reading the whole gene count. A module
        reads as the share of its families that lineage ``k`` carries.

        It needs no horizon breakpoint, and that is what makes the race exact rather than thinned:
        gene content changes only when a genome event fires, and an event ends the current step, so
        every rate is already constant between two events."""
        if fid is None:
            return len(gen[k])
        if isinstance(fid, tuple):
            return sum(counts.holds(k, f) for f in fid) / len(fid)
        return "present" if counts.holds(k, fid) else "absent"
    # the per-family weight sums, carried across events rather than rebuilt each time (see the class).
    # The families that wrote their own rate ride in the same structure under a suffixed key, because
    # summing a table over a lineage's copies is the same work whichever kind of number is in it.
    _tables = dict(fam_mult)
    if fam_fixed is not None:
        _tables.update({key + _FIXED: m for key, m in fam_fixed.items()})
    weights = _FamilyWeights(_tables, gen) if any_family else None
    counts = _FamilyCounts(gen)      # the family cap's question, answered without walking a genome

    # Recipient rules that read live gene content get their reader now, when the family ids, the
    # genomes and the copy counts all exist. A family's own rule is keyed by the id the family was
    # given, because that id is what a copy carries when it is picked to move.
    def live_reader(rule):
        return _LiveGeneContent(live_target(rule.driver, named, module_map or {}), gen, pos, counts)

    if isinstance(transfer_to, Driven) and names_a_live_level(transfer_to.driver):
        to_traj = live_reader(transfer_to)
    fam_choice: dict[int, tuple] = {}
    for i, (rule, groups_i, traj_i) in fam_choice_prepared.items():
        if isinstance(rule, Driven) and names_a_live_level(rule.driver):
            traj_i = live_reader(rule)
        fam_choice[named[declared[i].name]] = (rule, groups_i, traj_i)

    # four bare numbers — the per-copy trio and a per-lineage origination, no modifier on any rate,
    # no family writing its own — is the common run, and it needs none of the loop's context
    # machinery: each total is scope(base) exactly, resolved here once. A rate whose base is None
    # carries a set_by, which is a modifier, so `plain` is False and its 0.0 is never read.
    plain = (not (dup.modifiers or los.modifiers or org.modifiers or tra.modifiers)
             and fam_fixed is None and not any_per_lineage
             and dup.scope is PerCopy and los.scope is PerCopy and tra.scope is PerCopy
             and org.scope is PerLineage)
    dup_base, los_base = dup.base or 0.0, los.base or 0.0
    tra_base, org_base = tra.base or 0.0, org.base or 0.0

    # 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
        if live_keys and n > MAX_LIVE_COPIES:
            raise RuntimeError(
                f"the run passed {MAX_LIVE_COPIES} live gene copies at time {t:.3g} and is still "
                f"growing — a rate reading the genome's own content is feeding itself. Lower the "
                f"rates, flatten the mapping the driver is read through, or set a max_family_size.")
        k_alive = len(alive)
        can_xfer = n > 0 and (k_alive >= 2 or self_transfer)  # a recipient must be able to exist
        next_species = schedule[si][0]  # the tree's own next event: who is alive changes only here
        # a family placed by `origins=` originates at a fixed instant, so it joins the horizon like
        # any other breakpoint: the waiting time can never step over it
        next_start = starts[start_i][0] if start_i < len(starts) else math.inf
        w_dup = w_los = w_org = w_tra = None
        if plain:
            # no modifier on any rate: each total is scope(base) exactly — the per-copy trio times
            # the live copies, origination times the living lineages — none of them ever changes on
            # its own (next_change is inf), and nothing below reads a context or a weight.
            r_dup = dup_base * n
            r_los = los_base * n
            r_org = org_base * k_alive
            r_tra = tra_base * n if can_xfer else 0.0
            horizon = min(next_species, next_start)
        else:
            ctx = {"copies": n, "lineages": k_alive, "time": t}
            # A copy-consuming event counted *per lineage* is counted per lineage that HOLDS a copy:
            # an empty genome offers nothing to duplicate or lose, so it must not contribute its
            # share of the total and then be picked with no victim inside it. Origination keeps
            # `ctx` — an empty genome can still gain a family. Computed only when some rate needs
            # it, so the per-copy path does exactly the work it did before.
            if any_per_lineage:
                n_hosts = sum(1 for g in gen if g)
                host_ctx = {"copies": n, "lineages": n_hosts, "time": t}
            else:
                n_hosts, host_ctx = 0, ctx
            # A rate is per-lineage when something makes it differ between lineages, and two things
            # do: a driver, whose value is the one on that branch, and a per-lineage draw, whose
            # factor was drawn for that branch. Either way its effective rate is summed over the
            # living lineages (each read with its own copy count) and the weights are kept for the
            # affected-lineage pick — the species_tree._grow shape. A rate with neither stays pooled
            # (one .effective, uniform pick), so a run carrying neither is byte-identical to before.
            # For transfer the affected lineage is the donor, so both weight who donates. A run
            # carrying a per-family draw as well multiplies the two — the branch's factor is the
            # lineage's, the multipliers are its contents' — which is what `_lineage_weights` does
            # with `fam_sums`.
            fw = None
            if any_driven:  # each lineage's driver values, read before the weights that multiply them in
                drivers = [{**{key: trajs[key].value(alive[k], t) for key in trajs},
                            **{src: live_value(fid, k) for src, fid in live_reads}}
                           for k in range(k_alive)]
            if any_lineage:  # each lineage's drawn factor, in `alive` order, for the same weights
                lin_at = {key: [lin_mult[key][alive[k]] for k in range(k_alive)]
                          for key in lin_by if lin_by[key]}
            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.
                assert weights is not None       # `any_family` is exactly when it was built
                fw = weights.current(gen)
                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}
                own_sums = (lambda key: fw[key + _FIXED]) if fam_fixed is not None else (lambda key: None)
                if any_family_driven:
                    # A family whose own rate carries a verb has a different rate on each lineage, so its
                    # share is summed here at each event (its copies on lineage k times its rate read on
                    # k) and added to the fixed families' sums, which `_FamilyWeights` carries.
                    fixed_sums = own_sums
                    with_driven = {
                        key: [total + sum(counts.count(k, f) * family_rate_on(rate, k)
                                          for f, rate in fam_driven_ids[key].items()
                                          if counts.holds(k, f))
                              for k, total in enumerate(fw[key + _FIXED])]
                        for key in fam_driven}

                    def own_sums(key, _fixed=fixed_sums, _driven=with_driven):
                        return _driven[key] if key in _driven else _fixed(key)

                def unit_at(key, k, _rates={"duplication": dup, "loss": los, "transfer": tra}):
                    """The run's unit rate as lineage ``k`` reads it. Identical to ``unit[key]`` unless
                    the rate is driven or varies among lineages, and then it is the number
                    `_lineage_weights` used for that lineage — which the copy pick has to use too, or
                    the totals and the pick disagree about how a written family rate compares with
                    the run's."""
                    factor = lin_at[key][k] if lin_by[key] else 1.0
                    if not any_driven:
                        return unit[key] * factor
                    return _rates[key].effective(copies=1, lineages=1, time=t, drivers=drivers[k],
                                                 carried_factor=factor)
                w_dup = _family_weights(unit["duplication"], fw["duplication"],
                                        own_sums("duplication"),
                                        lin_at["duplication"] if lin_by["duplication"] else None)
                w_los = _family_weights(unit["loss"], fw["loss"], own_sums("loss"),
                                        lin_at["loss"] if lin_by["loss"] else None)
                if can_xfer:
                    w_tra = _family_weights(unit["transfer"], fw["transfer"], own_sums("transfer"),
                                            lin_at["transfer"] if lin_by["transfer"] else None)
            if any_driven or any_lineage:
                # each rate reads the driver values only if it is driven and the drawn factors only
                # if it carries a draw, so a rate with one of the two is unchanged by the other
                def weigh(key, rate, mods):
                    return _lineage_weights(
                        rate, gen, k_alive, t, drivers if mods else None,
                        fw[key] if fw is not None else None,
                        own_sums(key) if fw is not None else None,
                        lin_at[key] if lin_by[key] else None)
                if dup_mods or lin_by["duplication"]:
                    w_dup = weigh("duplication", dup, dup_mods)
                if los_mods or lin_by["loss"]:
                    w_los = weigh("loss", los, los_mods)
                # Origination is written out rather than taken through `_lineage_weights`, which
                # zeroes an empty genome: that is right for an event that consumes a copy and wrong
                # here, because a lineage holding nothing can still gain a family. It needs no
                # fam_sums branch either — origination can carry no per-family draw (refused above:
                # when it is read there is no family yet to have drawn one).
                if org_mods and lin_by["origination"]:
                    w_org = [org.effective(copies=len(gen[k]), lineages=1, time=t,
                                           drivers=drivers[k],
                                           carried_factor=lin_at["origination"][k])
                             for k in range(k_alive)]
                elif org_mods:
                    w_org = [org.effective(copies=len(gen[k]), lineages=1, time=t, drivers=drivers[k])
                             for k in range(k_alive)]
                elif lin_by["origination"]:
                    w_org = [org.effective(copies=len(gen[k]), lineages=1, time=t,
                                           carried_factor=lin_at["origination"][k])
                             for k in range(k_alive)]
                if (tra_mods or lin_by["transfer"]) and can_xfer:
                    w_tra = weigh("transfer", tra, tra_mods)
            r_dup = sum(w_dup) if w_dup is not None else (
                dup.effective(**(host_ctx if dup_per_lineage else ctx)) if n else 0.0)
            r_los = sum(w_los) if w_los is not None else (
                los.effective(**(host_ctx if los_per_lineage else 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(**(host_ctx if tra_per_lineage else ctx)) if can_xfer else 0.0)
            horizon = min(next_species, next_start, dup.next_change(t), los.next_change(t),
                          org.next_change(t), tra.next_change(t))
            if any_family_driven:  # a family's own changing_at moves its rate too, so step there
                horizon = min(horizon, min(rate.next_change(t) for table in fam_driven.values()
                                           for rate in table.values()))
            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)
        total = r_dup + r_los + r_org + r_tra

        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"],
                                                  own_rates("duplication", k),
                                                  unit_at("duplication", k))
                             if any_family else int(rng.integers(len(gen[k]))))
                    elif dup_per_lineage:  # every occupied genome equally likely, then a copy in it
                        k = _pick_host(rng, gen, n_hosts)
                        j = int(rng.integers(len(gen[k])))
                    else:
                        k, j = _pick_copy(rng, gen, n)
                    fam = gen[k][j].family
                    if not counts.at_cap(k, fam, cap):
                        _duplicate(gen[k], j, tree.nodes[alive[k]], t, events, new_copy)
                        counts.added(k, fam)
                        total_copies += 1
                        if weights is not None:
                            weights.touched(k)
                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"],
                                                  own_rates("loss", k), unit_at("loss", k))
                             if any_family else int(rng.integers(len(gen[k]))))
                    elif los_per_lineage:
                        k = _pick_host(rng, gen, n_hosts)
                        j = int(rng.integers(len(gen[k])))
                    else:
                        k, j = _pick_copy(rng, gen, n)
                    counts.removed(k, gen[k][j].family)      # before the copy leaves the genome
                    _lose_at(gen[k], j, tree.nodes[alive[k]], t, events)
                    total_copies -= 1
                    if weights is not None:
                        weights.touched(k)
                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)
                    counts.added(k, gen[k][-1].family)       # the copy _originate just appended
                    total_copies += 1
                    if weights is not None:
                        weights.touched(k)
                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"],
                                                   own_rates("transfer", kd),
                                                   unit_at("transfer", kd))
                              if any_family else int(rng.integers(len(gen[kd]))))
                    elif tra_per_lineage:  # every occupied genome donates equally often
                        kd = _pick_host(rng, gen, n_hosts)
                        jd = int(rng.integers(len(gen[kd])))
                    else:
                        kd, jd = _pick_copy(rng, gen, n)
                    choice, choice_groups, choice_traj = transfer_to, group_of, to_traj
                    if fam_choice:  # the moving copy's family may have a recipient rule of its own
                        own_choice = fam_choice.get(gen[kd][jd].family)
                        if own_choice is not None:
                            choice, choice_groups, choice_traj = own_choice
                    delta, kr = _do_transfer(rng, tree, alive, gen, counts, kd, jd, t, events,
                                             new_copy, choice, replacement, self_transfer,
                                             depth, choice_traj, cap, choice_groups)
                    total_copies += delta
                    if weights is not None and kr is not None:
                        weights.touched(kr)   # only the recipient's composition changed (see there)
                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]
                k_out = pos[i]
                g = gen[k_out]
                genomes[i] = tuple(g)  # finalise this lineage (extant, extinct, or unsampled)
                total_copies -= len(g)
                retire(alive, gen, pos, k_out)
                inherited = counts.retired(k_out)   # what the daughters below inherit, if any
                if weights is not None:
                    weights.retired(k_out)
                node = tree.nodes[i]
                if node.children:  # a speciation: each gene re-ids into each daughter
                    per_daughter = []
                    for c in node.children:
                        child_genome, rows = [], []
                        for old in g:  # ZOMBI1: the gene ends here and continues under a fresh id
                            nc = new_copy(old.family)
                            child_genome.append(nc)
                            rows.append(GeneEdge(t, "speciation", c, old.family, nc.id, parent=old.id))
                        per_daughter.append(rows)
                        enter(alive, gen, pos, c, child_genome)
                        counts.entered_like(inherited)   # a re-id of the parent: same families
                        if weights is not None:
                            weights.entered(child_genome)
                        total_copies += len(child_genome)
                    # the ids are minted daughter by daughter (which is what fixes them), but a gene's
                    # two rows are recorded together: one gene ending is one event, and the log writes
                    # it as one row with both daughters in it.
                    for pair in zip(*per_daughter):
                        events.extend(pair)
                si += 1
        elif start_i < len(starts) and horizon == next_start:
            # a placed family arrives. The lineage is live by construction (its time was checked
            # against that branch's own life), and a tie with the tree's schedule falls to the
            # branch above: the daughters have entered by the time this runs.
            t = horizon
            while start_i < len(starts) and starts[start_i][0] == t:
                _, lineage, fam = starts[start_i]
                k = pos[lineage]
                c = new_copy(fam)
                gen[k].append(c)
                events.append(GeneEdge(t, "origination", lineage, fam, c.id))
                counts.added(k, fam)
                total_copies += 1
                if weights is not None:
                    weights.touched(k)
                start_i += 1
        else:
            t = horizon  # a skyline breakpoint: advance and re-evaluate the (now changed) rate

    bar.close()
    links = links_of({"duplication": dup, "transfer": tra, "loss": los, "origination": org},
                     transfer_to, declared, fam_driven, fam_transfer_to, module_map or {})
    multipliers = (multipliers_of(fam_mult, FAMILY_TARGETS, own_cells)
                   if any(fam_by.values()) else {})
    per_lineage_multipliers = (lineage_multipliers_of(lin_mult, tree.labels())
                               if any_lineage else {})
    return FamilyGenomesResult(tree, genomes, events, seed, named, module_map, initial_genome,
                               cap, links, multipliers, per_lineage_multipliers)

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, families=None, joint=False, max_family_size=10, seed=None, stream_to=None, outputs=None, compact=False, bundle=None, progress=False, **retired) -> 'OrderedGenomesResult | StreamedRun'

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 events.

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 extent 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 rates fission/fusion/chromosome_loss are 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; families=[family("toxin")] additionally declares named families (remembered in result.family_names for result.has_family(node, "toxin")), as in the family core; replacement / self_transfer behave as in the family core. So does transfer_to, which chooses who receives — "uniform", "distance" / Distance(decay=) (closer relatives likelier), Clades({...}, Between({...})) (weight by the donor's and recipient's named clade) or Recipients().weighted_by(driver, mapping) (weight by another level; see below). What moves is a block of genes rather than a single copy, and the block arrives whole, so the rule chooses the recipient lineage exactly as it does at the family resolution.

The chromosome events change chromosome number: fission (split), fusion (merge, between two chromosomes of the same topology — the reticulation; a ring and a molecule with two ends cannot become one molecule, so a genome of one of each never fuses), 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.

A family placed by hand. origins=[("n5", 0.4)] originates a family on lineage n5 at time 0.4 — the ordinary origination event, at a point you choose rather than one that is drawn, and here too the founding gene lands on a uniformly-chosen chromosome and position. It adds to whatever initial_families and origination already give you. See resolve_origins.

Conditioning (a trait drives a rate). Any rate here may be driven by another level — inversion = PerCopy(0.3).scaled_by(habitat, {"host": 4.0, "free": 1.0}) scales each lineage's inversion rate by the habitat on that branch, read from a trait grown first (the finished TraitsResult, or the trait_events.tsv it wrote). A driven rate is then per lineage: it is summed over the living lineages, each read with its own gene count, chromosome count and driver value; the lineage an event lands on is drawn with those same weights, and the gene inside it uniformly, because the gene count is already in the weight. The Gillespie steps at every mid-branch switch of the driver rather than averaging over a branch (SPEC §2). For transfer the driven lineage is the donor, so a driven transfer says how often a lineage donates.

Conditioning (a trait drives who receives). transfer_to = Recipients().weighted_by(driver, 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, a weight, not a rate). Candidate lineage k gets weight mapping(driver value on k now) and receives with probability w_k / Σw. Weight 0 means "cannot receive"; when every candidate weighs 0 the transfer does not happen at all, and the donor's chromosome is left untouched. A Between({...}) mapping reads the donor's value too, so transfer can be steered between guilds; Clades({...}, Between({...})) is the same steering by named clade, read off the tree instead of a driver. Because a weight is not a rate, a driven transfer_to adds no Gillespie breakpoint and composes freely with a driven transfer rate.

Conditioning (a trait drives an extent). An extent takes the same modifiers a rate does (SPEC §6) — inversion_extent = Extent(4).scaled_by(habitat, {"host": 3.0, "free": 1.0}) makes a host-restricted lineage invert longer runs of genes, which is a different statement from raising its inversion rate. An extent's modifier is read at the instant an event fires, so it changes how much a run takes and never how often one starts, and it adds no Gillespie breakpoint.

Joint runs (gene content drives a rate). With joint=True a rate, an extent or transfer_to can read the gene content the run is building: "genomes:<family>" (whether a declared family is there), "genomes:count" (how many genes the lineage has) or "genomes:module:<group>" (the share of a module's families it carries), read on each lineage as the run goes. On a rate or an extent the factor belongs to the lineage, so it composes with extents as a trait's does (SPEC §6). On transfer_to each candidate is weighted by its own genome when the transfer fires.

A family's own rate. families=[family("B", loss=0.8)] gives one family its own duplication, transfer or loss, fixed or read from a driver, as at the family resolution. It applies to the segment, as a per-family draw does (SPEC §6): every gene carries its family's rate or the run's, an event's total is their sum, and a run of genes is chosen by the mean rate of the genes it covers. With extents longer than one gene, a family's own rate decides where events start, not which genes they take: a family whose rate is 0 can still be removed by a segment that starts on a neighbour.

A family's own transfer_to. family("B", transfer_to=Recipients().weighted_by("genomes:A", {"present": 20.0, "absent": 1.0})) gives one family its own recipient rule, in any form the run's transfer_to takes, as at the family resolution. A transferred segment arrives whole, and every gene it covers carries its family's rule or the run's. A lineage can receive the segment only where every one of those rules allows: their weights multiply, each rule once, before they are normalised. With the run's rule uniform, a segment carrying B keeps B's preference exactly, unless it also carries a family with a rule of its own. A weight of 0 in any rule excludes that lineage, and when no lineage is left the transfer does not happen.

a per-family draw and a driven rate cannot be set in the same run: one weights lineages by a driver and the other weights the segment by what it covers.

Streaming to disk. stream_to=DIR writes the run's files as the run goes and keeps none of its record in memory, for a run whose record would not fit. The run is the same run: the same seed gives the same events, and the files hold what result.write(DIR) writes for it. Two things differ. gene_order.tsv lists each node's rows when its branch ends rather than in node order, and a StreamedRun comes back instead of an OrderedGenomesResult. The gene trees are built at the end from the event log, one group of families at a time. outputs= picks the files, as write() takes them, and without stream_to it is an error. compact=True writes the gene trees as one gene_trees.tsv, as write(compact=True) does, and is an error without stream_to too. The genomes still alive are kept, as the run needs them, so at the end every extant genome is in memory at once.

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, families=None, joint=False,
                             max_family_size=10, seed=None, stream_to=None, outputs=None,
                             compact=False, bundle=None, progress=False,
                             **retired) -> "OrderedGenomesResult | StreamedRun":
    """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 events.

    **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 **extent** 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 rates ``fission``/``fusion``/``chromosome_loss``
    are **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**; ``families=[family("toxin")]`` additionally
    declares **named** families (remembered in ``result.family_names`` for ``result.has_family(node,
    "toxin")``), as in the family core; ``replacement`` / ``self_transfer`` behave as in the family
    core. So does ``transfer_to``, which **chooses who receives** — ``"uniform"``,
    ``"distance"`` / ``Distance(decay=)`` (closer relatives likelier), ``Clades({...}, Between({...}))``
    (weight by the donor's and recipient's named clade) or ``Recipients().weighted_by(driver, mapping)`` (weight by
    another level; see below). What moves is a block of genes rather than a single copy, and the block
    arrives whole, so the rule chooses the recipient lineage exactly as it does at the family
    resolution.

    The **chromosome events** change chromosome *number*: ``fission`` (split), ``fusion`` (merge,
    between two chromosomes of the **same topology** — the reticulation; a ring and a molecule with
    two ends cannot become one molecule, so a genome of one of each never fuses),
    ``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``.

    **A family placed by hand.** ``origins=[("n5", 0.4)]`` originates a family on lineage ``n5`` at
    time ``0.4`` — the ordinary origination event, at a point you choose rather than one that is
    drawn, and here too the founding gene lands on a uniformly-chosen chromosome and position. It
    adds to whatever ``initial_families`` and ``origination`` already give you. See
    `resolve_origins`.

    **Conditioning (a trait drives a rate).** Any rate here may be *driven by another level* —
    ``inversion = PerCopy(0.3).scaled_by(habitat, {"host": 4.0, "free": 1.0})`` scales each lineage's
    inversion rate by the habitat on that branch, read from a trait grown first (the finished
    ``TraitsResult``, or the ``trait_events.tsv`` it wrote). A driven rate is then *per lineage*: it is
    summed over the living lineages, each read with its own gene count, chromosome count and driver
    value; the lineage an event lands on is drawn with those same weights, and the gene inside it
    uniformly, because the gene count is already in the weight. The Gillespie steps at **every**
    mid-branch switch of the driver rather than averaging over a branch (SPEC §2). For ``transfer``
    the driven lineage is the **donor**, so a driven ``transfer`` says how often a lineage *donates*.

    **Conditioning (a trait drives who receives).** ``transfer_to =
    Recipients().weighted_by(driver, 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, a weight, not a rate). Candidate lineage ``k`` gets weight ``mapping(driver value on k now)``
    and receives with probability ``w_k / Σw``. Weight 0 means "cannot receive"; when every candidate
    weighs 0 the transfer does not happen at all, and the donor's chromosome is left untouched. A
    ``Between({...})`` mapping reads the **donor's** value too, so transfer can be steered between
    guilds; ``Clades({...}, Between({...}))`` is the same steering by named clade, read off the tree
    instead of a driver. Because a weight is not a rate, a driven ``transfer_to`` adds no Gillespie
    breakpoint and composes freely with a driven ``transfer`` rate.

    **Conditioning (a trait drives an extent).** An extent takes the same modifiers a rate does
    (SPEC §6) — ``inversion_extent = Extent(4).scaled_by(habitat, {"host": 3.0, "free": 1.0})`` makes a
    host-restricted lineage invert *longer runs of genes*, which is a different statement from raising
    its inversion rate. An extent's modifier is read at the instant an event fires, so it changes how
    much a run takes and never how often one starts, and it adds no Gillespie breakpoint.

    **Joint runs (gene content drives a rate).** With ``joint=True`` a rate, an extent or
    ``transfer_to`` can read the gene content the run is building: ``"genomes:<family>"`` (whether a
    declared family is there), ``"genomes:count"`` (how many genes the lineage has) or
    ``"genomes:module:<group>"`` (the share of a module's families it carries), read on each lineage
    as the run goes. On a rate or an extent the factor belongs to the lineage, so it composes with
    extents as a trait's does (SPEC §6). On ``transfer_to`` each candidate is weighted by its own
    genome when the transfer fires.

    **A family's own rate.** ``families=[family("B", loss=0.8)]`` gives one family its own
    ``duplication``, ``transfer`` or ``loss``, fixed or read from a driver, as at the family
    resolution. It applies to the segment, as a per-family draw does (SPEC §6): every gene carries its
    family's rate or the run's, an event's total is their sum, and a run of genes is chosen by the mean
    rate of the genes it covers. With extents longer than one gene, a family's own rate decides where
    events start, not which genes they take: a family whose rate is 0 can still be removed by a
    segment that starts on a neighbour.

    **A family's own transfer_to.** ``family("B", transfer_to=Recipients().weighted_by("genomes:A",
    {"present": 20.0, "absent": 1.0}))`` gives one family its own recipient rule, in any form the
    run's ``transfer_to`` takes, as at the family resolution. A transferred segment arrives whole, and
    every gene it covers carries its family's rule or the run's. A lineage can receive the segment
    only where every one of those rules allows: their weights multiply, each rule once, before they
    are normalised. With the run's rule uniform, a segment carrying ``B`` keeps ``B``'s preference
    exactly, unless it also carries a family with a rule of its own. A weight of 0 in any rule
    excludes that lineage, and when no lineage is left the transfer does not happen.

    a per-family draw and a driven rate cannot be set in the same run: one weights lineages by a driver and
    the other weights the segment by what it covers.

    **Streaming to disk.** ``stream_to=DIR`` writes the run's files as the run goes and keeps none of
    its record in memory, for a run whose record would not fit. The run is the same run: the same seed
    gives the same events, and the files hold what ``result.write(DIR)`` writes for it. Two things
    differ. ``gene_order.tsv`` lists each node's rows when its branch ends rather than in node order,
    and a `StreamedRun` comes back instead of an `OrderedGenomesResult`. The gene trees are built at the
    end from the event log, one group of families at a time. ``outputs=`` picks the files, as
    ``write()`` takes them, and without ``stream_to`` it is an error. ``compact=True`` writes the
    gene trees as one ``gene_trees.tsv``, as ``write(compact=True)`` does, and is an error without
    ``stream_to`` too. The genomes still alive are kept, as the run needs them, so at the end every
    extant genome is in memory at once.
    """
    compact = compact_or_bundle(compact, bundle, "simulate_genomes_ordered")   # its 0.52.0 name
    tree = as_tree(tree, level="genomes")
    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 compact and stream_to is None:
        raise ValueError(
            "compact applies to a streamed run (stream_to=DIR), which writes the files itself; for "
            "an in-memory run pass it when you call result.write(compact=True).")
    if stream_to is not None:
        outputs = tuple(OrderedGenomesResult.OUTPUTS if outputs is None else outputs)
        if unknown := [o for o in outputs if o not in OrderedGenomesResult.OUTPUTS]:
            raise ValueError(f"unknown stream outputs {unknown}; choose from "
                             f"{list(OrderedGenomesResult.OUTPUTS)}")
    labels = _topologies(chromosomes, topology)
    n_initial_chrom = chromosomes
    # this slice implements each event's default scope and the cells IMPLEMENTED_MODIFIERS
    # declares: changing_at (skyline), scaled_by (a conditioned/joint driver, per lineage), set_by (a
    # driver that replaces the base), a per-family draw and a draw among lineages —
    # the per-family one with the weight on the SEGMENT rather than on its starting gene (SPEC §6,
    # and _pick_run_by_family), the per-lineage one on the lineage, as a driver's factor is. The
    # per-lineage draw reaches every rate here, the chromosome events included: it scales whatever
    # the event acts on, where a per-family weight has to reach the genes a segment covers.
    _rates: dict[str, Rate] = {}
    for label, spec, want in (("duplication", duplication, PerCopy), ("transfer", transfer, PerCopy),
                              ("loss", loss, PerCopy), ("origination", origination, PerLineage),
                              ("inversion", inversion, PerCopy),
                              ("transposition", transposition, PerCopy),
                              ("translocation", translocation, PerCopy),
                              ("fission", fission, PerChromosome), ("fusion", fusion, PerChromosome),
                              ("chromosome_origination", chromosome_origination, PerLineage),
                              ("chromosome_loss", chromosome_loss, PerChromosome)):
        rate = as_rate(spec, default_scope=want)
        # An event that acts on **genes** takes either answer to *per what?*: per copy (the default —
        # each gene independently at risk, so a bigger genome turns over faster) or per lineage (a
        # fixed budget, the same however much the genome holds). The chromosome rates and origination
        # keep one scope each: origination creates families, so per copy it would be base × 0 in an
        # empty genome, and the chromosome rates are not implemented per lineage.
        legal = (want, PerLineage) if want is PerCopy else (want,)
        # `rate.scope` holds the scope CLASS, not an instance — a scope constructor returns the rate
        # itself — so this is an identity test against the legal set rather than an isinstance one.
        assert rate.scope is not None            # as_rate fills the level's default where none was written
        if rate.scope not in legal:
            raise ValueError(
                f"{label} has a {rate.scope.__name__} scope, but the ordered genome engine "
                f"takes {' or '.join(s.__name__ for s in legal)} for {label}."
            )
        for m in rate.modifiers:
            if m.reads == (DRAWN, "families") and label == "origination":
                raise ValueError(
                    "origination carries a per-family draw, 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. "
                    "Write varying_among('families', …) on duplication, transfer, loss, inversion, "
                    "transposition or translocation; writing one such object on several of them gives "
                    "a family-wide tempo, since one object is one draw.")
            if m.reads == (DRAWN, "families") and rate.scope is not PerCopy:
                raise ValueError(
                    f"{label} carries a per-family draw on a {rate.scope.__name__} scope. A per-family "
                    f"weight has to reach the genes an event covers, so it applies to the per-copy "
                    f"gene events only — not to the chromosome events, which act on whole replicons.")
            if isinstance(m, Driven):
                check_not_a_kernel(m.mapping, label=label)
            if not is_implemented(m, IMPLEMENTED_MODIFIERS, "genomes.ordered"):
                raise ValueError(
                    f"{label} carries {describe(m)}, which the ordered genome engine does not "
                    f"support. It takes changing_at (skyline), scaled_by (a conditioned or joint "
                    f"driver), set_by (a driver that replaces the base), varying_among('families', "
                    f"…) (per-family heterogeneity, weighted on the segment an event covers) and "
                    f"varying_among('lineages', …) (one multiplier per species branch, drawn or "
                    f"inherited)."
                )
        _rates[label] = rate
    # the eleven rates keep short names in the Gillespie loop below; the dict is what the driver
    # resolution and the per-lineage weights walk, so neither has to name all eleven again
    dup, tra, los, org = (_rates["duplication"], _rates["transfer"], _rates["loss"],
                          _rates["origination"])
    inv, trp, trl = _rates["inversion"], _rates["transposition"], _rates["translocation"]
    fis, fus = _rates["fission"], _rates["fusion"]
    cor, clo = _rates["chromosome_origination"], _rates["chromosome_loss"]
    for label, r in _rates.items():
        r.check_one_base(label)
        # 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 branch's multiplier to be independent
        # of its parent's and inherited from it at once — there is no model there to implement.
        check_one_memory(tuple(m for m, _ in r.carried_modifiers(unit="lineages")),
                         label=label, unit="lineages")
    # Over the whole RUN, not per rate, and getting that wrong was a real bug: a per-family draw
    # anywhere makes the engine take its per-family path for **every** gene rate, summing each one
    # over the live genes — so a `PerLineage` rate elsewhere in the same run had its total counted
    # per copy while its acting lineage was still drawn uniformly among occupied genomes. The total
    # and the pick then said different things, which is the one failure this engine must not have.
    _GENE_EVENTS = ("duplication", "transfer", "loss", "inversion", "transposition", "translocation")
    per_lineage_here = [lbl for lbl in _GENE_EVENTS if _rates[lbl].scope is PerLineage]
    drawn_here = [lbl for lbl in _GENE_EVENTS
                  if any(m.reads == (DRAWN, "families") for m in _rates[lbl].modifiers)]
    if per_lineage_here and drawn_here:
        raise ValueError(
            f"{', '.join(per_lineage_here)} is PerLineage while {', '.join(drawn_here)} carries a "
            f"per-family draw, and the two cannot share a run. Under PerCopy a family's multiplier "
            f"scales each gene's rate, so it changes the lineage's total; under PerLineage the total "
            f"is fixed whatever the genome holds, so the multiplier could only choose which segment "
            f"is taken. Those are different models and the choice is not made yet — write PerCopy "
            f"throughout for the first, or drop the per-family draw for the second.")
    if any(m.reads == (DRAWN, "families") for r in _rates.values() for m in r.modifiers) and \
            any(isinstance(m, Driven) for r in _rates.values() for m in r.modifiers):
        raise ValueError(
            "a per-family draw and a driver on the same run is not wired at the ordered resolution: "
            "a driver weights the lineage, and here a per-family draw has to weight the SEGMENT an "
            "event covers rather than the gene it started from, so the two are not one "
            "multiplication. The family resolution runs the pair — there a family's multiplier is "
            "the copy's, and the weight is simply the product. Use it, or use one of the two here.")
    # 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 a rate takes at this resolution minus a per-family draw (see
        `IMPLEMENTED_EXTENT_MODIFIERS`), and they scale the size drawn — ``changing_at`` in time,
        ``scaled_by`` on the lineage the event lands on."""
        e = as_extent(spec)
        rate_slot = label.removesuffix("_extent")
        for m in e.modifiers:
            if isinstance(m, Driven):
                check_not_a_kernel(m.mapping, label=label)
            if m.reads == (DRAWN, "families"):
                raise ValueError(
                    f"{label} carries a per-family draw, which an extent cannot mean: the size is drawn before "
                    f"the run's genes are known, and a run covers several families, so there is no "
                    f"one family to draw a factor for. Put it on {rate_slot}, where it weights "
                    f"the segment by what it covers.")
            if not is_implemented(m, IMPLEMENTED_EXTENT_MODIFIERS, "genomes.ordered"):
                raise ValueError(
                    f"{label} carries {describe(m)}, which the ordered genome engine does not "
                    f"support on an extent — it takes "
                    f"{', '.join(cell_name(w) for w in IMPLEMENTED_EXTENT_MODIFIERS)}.")
        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"))
    _extents: dict[str, Extent] = {
        "duplication_extent": dup_ext, "loss_extent": los_ext, "transfer_extent": tra_ext,
        "inversion_extent": inv_ext, "transposition_extent": trp_ext,
        "translocation_extent": trl_ext}
    if not 0.0 <= inversion_probability <= 1.0:
        raise ValueError(f"inversion_probability must be in [0, 1], got {inversion_probability!r}")
    # the choice (SPEC §5), validated in the one place all three resolutions share: the mapping's
    # numbers are weights over the candidate recipients, never a rate multiplier
    transfer_to = resolve_transfer_to(transfer_to)
    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}")
    check_no_retired_keywords(retired, where="simulate_genomes_ordered")
    # the family resolution's own resolver, so the two engines cannot disagree about what a
    # declaration means
    declared, module_map, started_named = resolve_families(families, tree)
    family_names = [f.name for f in declared]
    # A family's own duplication, transfer or loss, read by the family resolution's resolver. Here it
    # applies to the segment, as a per-family draw does (SPEC §6): every gene carries its family's rate
    # or the run's, and the extent stays the run's.
    fam_own, fam_driven_rates = resolve_family_rates(
        declared, {"duplication": dup, "transfer": tra, "loss": los})
    own_keys = sorted(set(fam_own) | set(fam_driven_rates))
    any_written = bool(own_keys)
    for key in own_keys:
        run_scope = _rates[key].scope
        assert run_scope is not None             # as_rate fills the level's default where none was written
        if run_scope is not PerCopy:
            raise ValueError(
                f"a family writes its own {key}, but the run's {key} is "
                f"{run_scope.__name__}. The two are summed over the same genes, so both are "
                f"counted per copy: write PerCopy for the run's {key}, or drop the family's.")
    # a family's own recipient rule, checked by the same resolver as the run's transfer_to
    fam_transfer_to = resolve_family_transfer_to(declared)

    # 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 a per-family draw 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)

    # Conditioning: a rate written with scaled_by reads a driver **per lineage**, so its rate stops being
    # one number for the whole live set and becomes one per lineage. Same machinery as the other two
    # resolutions — each driver 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 and no driven extent this is
    # empty and the loop stays exactly the pooled one, so an undriven run is untouched.
    driven = {label: [m for m in r.modifiers if isinstance(m, Driven)]
              for label, r in _rates.items()}
    ext_driven = {label: [m for m in e.modifiers if isinstance(m, Driven)]
                  for label, e in _extents.items()}
    # A **live** driver names gene content this run is building ("genomes:<family>", "genomes:count",
    # "genomes:module:<group>"), which makes the run joint. It is read off the genomes as the loop goes
    # rather than resolved into a trajectory, and its names are checked by the family resolution's own
    # resolver, so the two engines accept and refuse the same spellings.
    # a family's own rate reads its drivers as a run rate does, live or finished
    fam_rate_driven = [m for table in fam_driven_rates.values() for rate in table.values()
                       for m in rate.modifiers if isinstance(m, Driven)]
    live_rate_mods = [m for mods in (*driven.values(), fam_rate_driven) for m in mods
                      if names_a_live_level(m.driver)]
    live_ext_mods = [m for mods in ext_driven.values() for m in mods if names_a_live_level(m.driver)]
    # A recipient rule reading live gene content is read when a transfer fires. It is checked with the
    # rates' live drivers and makes the run joint too, but it moves no rate (see `prepare_transfer_to`).
    live_choices = [r for r in (transfer_to, *fam_transfer_to.values())
                    if isinstance(r, Driven) and names_a_live_level(r.driver)]
    resolve_live_drivers([*live_rate_mods, *live_ext_mods], set(family_names), joint=joint,
                         choice_mods=live_choices,
                         modules={name: len(members) for name, members in (module_map or {}).items()})
    by_key: dict = {}                   # driver key → its Driven (deduped: one driver resolves once)
    for mods in (*driven.values(), *ext_driven.values(), fam_rate_driven):
        for m in mods:
            if not names_a_live_level(m.driver):
                by_key.setdefault(m.key, m)
    resolved: dict = {}
    if by_key:
        resolved = {key: resolve_driver(m.driver, tree, step=m.step, path=m.path, level="genomes.ordered")
                    for key, m in by_key.items()}
        # a mapping whose states never occur leaves every lineage on the default factor, so the run
        # would secretly be the undriven model — refuse it here, naming the driver
        for mods in (*driven.values(), *ext_driven.values(), fam_rate_driven):
            for m in mods:
                if names_a_live_level(m.driver):
                    continue                    # checked above, and read live rather than resolved
                src = m.driver if isinstance(m.driver, str) else f"<{type(m.driver).__name__}>"
                check_mapping_fires(m.mapping, resolved[m.key].states(), driver_label=src)
    # 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(), fam_rate_driven) for m in mods}
    trajs = {key: traj for key, traj in resolved.items() if key in _rate_keys}
    # a live driver is a name by construction (names_a_live_level only admits strings)
    live_rate_keys = list(dict.fromkeys(cast(str, m.driver) for m in live_rate_mods))
    live_ext_keys = list(dict.fromkeys(cast(str, m.driver) for m in live_ext_mods))
    any_driven = bool(trajs) or bool(live_rate_keys)
    any_ext_driven = any(ext_driven.values())
    # The transfer_to slot is prepared **after** `trajs` is fixed, for the same reason: a driven
    # transfer_to is a weight, not a rate, so its trajectory must not join `trajs` and start adding
    # horizon breakpoints. `resolved` doubles as the driver cache, so a trait that drives both a rate
    # and who receives is loaded once and read from one trajectory.
    def prepare_choice(rule):
        """``(groups, trajectory)`` for a recipient rule. A rule reading live gene content has nothing
        to prepare here: its reader is built once the family ids and the genomes exist (below)."""
        if isinstance(rule, Driven) and names_a_live_level(rule.driver):
            return None, None
        return prepare_transfer_to(tree, rule, resolved, level="genomes.ordered")

    group_of, to_traj = prepare_choice(transfer_to)
    # A family's rule that equals the run's rule, or an earlier family's, is replaced by that rule, so
    # a segment carrying both counts it once (see `_do_transfer`). A name or a `Distance` is the same
    # rule when it is equal; any other rule only when it is the same object, as elsewhere in the
    # grammar one object read twice is one reading.
    choices: list[tuple] = [(transfer_to, group_of, to_traj)]
    fam_choice_prepared: dict[int, tuple] = {}
    for i, rule in fam_transfer_to.items():
        same = next((c for c in choices if c[0] is rule
                     or (isinstance(rule, (str, Distance)) and type(c[0]) is type(rule) and c[0] == rule)),
                    None)
        if same is None:
            same = (rule, *prepare_choice(rule))
            choices.append(same)
        fam_choice_prepared[i] = same

    rng, seed = stream("genomes", seed)     # own stream, and a drawn seed if none was given

    # Per-lineage multipliers: one per species branch, shared by every family that passes through it
    # (`zombi2.genomes.multipliers`). Drawn before any family is minted, because they come from the
    # tree, which exists before the run does; a run carrying none draws nothing, so it is
    # bit-identical to one from before this existed. On transfer the branch is the DONOR's, as a
    # driven transfer's is.
    lin_by = {label: tuple(m for m, _ in r.carried_modifiers(unit="lineages"))
              for label, r in _rates.items()}
    lin_mult = draw_lineage_multipliers(lin_by, tree, rng, targets=ORDERED_LINEAGE_TARGETS)
    any_lineage = bool(lin_mult)
    varying = {label: bool(mods) for label, mods in lin_by.items()}

    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 created and fixed for its whole life,
    # exactly as at the family resolution: one `Random` object read by two rates is one draw for
    # both, two objects are two draws. 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 some rate carries one.
    fam_by = {"duplication": tuple(m for m, _ in dup.carried_modifiers(unit="families")),
              "transfer": tuple(m for m, _ in tra.carried_modifiers(unit="families")),
              "loss": tuple(m for m, _ in los.carried_modifiers(unit="families")),
              "inversion": tuple(m for m, _ in inv.carried_modifiers(unit="families")),
              "transposition": tuple(m for m, _ in trp.carried_modifiers(unit="families")),
              "translocation": tuple(m for m, _ in trl.carried_modifiers(unit="families"))}
    any_family = any(fam_by.values())
    fam_mult: dict[str, _FamilyWeights] = {key: _FamilyWeights() for key in fam_by}

    # Which gene-level rates are a fixed per-lineage budget rather than a per-gene risk. Read once:
    # it decides both how the total is counted and how the acting lineage is picked, and those two
    # must never disagree.
    per_lineage = {label: _rates[label].scope is PerLineage
                   for label in ("duplication", "transfer", "loss",
                                 "inversion", "transposition", "translocation")}
    any_per_lineage = any(per_lineage.values())

    def new_family() -> int:
        nonlocal family_counter
        f = family_counter
        family_counter += 1
        if any_family:
            # one draw per distinct modifier object for this family, shared across its rates (see
            # `values_at_birth`): one object written on two rates is one number.
            shared: dict[int, float] = {}
            for key, mods in fam_by.items():
                fam_mult[key][f] = math.prod(values_at_birth(mods, rng, shared))
        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, ...]] = {}
    # A streamed run writes these four as it goes and keeps none of them (`_OrderedStream`): the engine
    # appends to the stream's own lists, and the stream writes and empties them in batches.
    to_disk = _OrderedStream(stream_to, outputs, tree, compact) if stream_to is not None else None
    events: list[GeneEdge] = to_disk.edges if to_disk is not None else []
    event_positions: list[EventPosition] = to_disk.positions if to_disk is not None else []
    rearrangements: list[Inversion | Transposition | Translocation] = (
        to_disk.rearrangements if to_disk is not None else [])
    chromosome_events: list[ChromosomeEvent] = (
        to_disk.chromosome_events if to_disk is not None else [])

    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, []))
        # `initial`, not `origination`: a replicon the run *starts* with is not something it did, so
        # counting `origination` in the log gives the de-novo replicons alone
        chromosome_events.append(ChromosomeEvent(t, "initial", root.id, (), (cid,)))
    # the origin'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]
        _live(chrom).append(new_gene(fam, +1))
        events.append(GeneEdge(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))
    # a family's own rate by the id the family is given: fixed ones as numbers, driven ones as the Rate
    fam_fixed_by_id: dict[str, dict[int, float]] = {key: {} for key in own_keys}
    fam_driven_by_id: dict[str, dict[int, _GeneRate]] = {key: {} for key in own_keys}
    named: dict[str, int] = {}  # a minted id per declared name, dealt round-robin after the anonymous ones
    named_starts: list[tuple[float, int, int]] = []
    for j, name in enumerate(family_names):
        fam = new_family()
        named[name] = fam
        for key in own_keys:
            if j in fam_own.get(key, {}):
                fam_fixed_by_id[key][fam] = fam_own[key][j]
            elif j in fam_driven_rates.get(key, {}):
                fam_driven_by_id[key][fam] = _GeneRate(fam_driven_rates[key][j])
        if j in started_named:
            # given an `origin`, so it arrives there rather than at the tree's origin — the same
            # event, at a point chosen instead of drawn
            t_p, lineage = started_named[j]
            named_starts.append((t_p, lineage, fam))
            continue
        chrom = initial_chroms[(initial_families + j) % n_initial_chrom]
        _live(chrom).append(new_gene(fam, +1))
        events.append(GeneEdge(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
    # the ids of the families `origins=` places: minted here, straight after the initial and named
    # ones and in the order they were written, so the same origins name the same families at either
    # resolution. Each starts at its own time, in the loop below.
    starts = sorted(named_starts)
    start_i = 0
    initial_genome = tuple(Chromosome(c.id, c.topology, list(c.genes)) for c in initial_chroms)
    enter(alive, gen, pos, root.id, initial_chroms)
    counts = _GeneCounts(gen)       # genes per family on every living lineage, changed with the genomes
    # What a step reads per lineage, kept across steps so a step rebuilds only the lineages an event
    # changed (`_LineageRows`). A run with no driven rate, no draw of either kind and no own rate
    # reads nothing per lineage, and then the rows are there only to count the living lineages.
    # A rate is read lineage by lineage when something makes it differ between lineages, and two
    # things do: a driver, whose value is the one on that branch, and a per-lineage draw, whose
    # multiplier was drawn for that branch. Either way the rate gets a `w` row, which the step totals
    # and an event draws the acting lineage by.
    lineage_rates = {label: rate for label, rate in _rates.items()
                     if driven[label] or varying[label]}
    # Where each `w` row comes from. A rate carrying a per-family draw as well keeps its weight in
    # `fw` — the family multipliers summed over its genes, which the segment pick needs unscaled —
    # so its row is that sum times the branch's multiplier, built beside `fw` rather than from the
    # rate alone. A driver and a per-family draw are refused in one run, so these two never overlap.
    w_from_family = {label for label in lineage_rates if any_family and fam_by.get(label)}
    w_from_rate = {label for label in lineage_rates if label not in w_from_family}
    any_rows = any_driven or any_family or any_written or any_lineage
    rows = _LineageRows(lineage_rates if (any_driven or any_lineage) else (),
                        fam_mult if any_family else (), own_keys if any_written else ())
    for genome in gen:
        rows.entered(genome)
    # Whether anything in this run changes with time on its own. A row reads time only through a
    # schedule — on a rate, on a family's own rate, or inside a driver's mapping — and through a
    # driver grown before the run. With none of those, the tree's own schedule moves `t` without
    # moving a single number a row holds, so a speciation marks the daughters that enter and leaves
    # every other lineage's row alone. That is the difference between a step at a speciation costing
    # the living lineages and costing two.
    # The rates that change on their own at some instant: a schedule on the rate, or on an entry of a
    # driver's mapping. `next_change` names the first change strictly after an instant, so a rate with
    # none after -inf has none at all. The loop sets its horizon from these alone, instead of asking
    # every rate at every step — with 200 declared families that was 211 questions a step, each
    # answered "never".
    timed_rates = [r for r in (dup, los, org, tra, inv, trp, trl, fis, fus, cor, clo,
                               *(fr for table in fam_driven_rates.values() for fr in table.values()))
                   if math.isfinite(r.next_change(-math.inf))]
    time_varying = bool(trajs) or bool(timed_rates)
    # one entry per event class whose families write their own rate, reading the rows as they stand:
    # the lineage's summed rates to draw it by, and `table(k)` for the rate of each of its genes
    own_pick = {key: (rows.own_sums[key],
                      _own_table(rows.own_units[key], rows.own_owned[key],
                                 fam_mult[key] if any_family else None, rows.own_largest[key]))
                for key in (own_keys if any_written else ())}
    # the run's own rate for each of those event classes, as the genes without an own rate carry it
    run_gene_rates = {key: _GeneRate(_rates[key]) for key in (own_keys if any_written else ())}
    # a family given an `origin` is not in the root genome — it arrives later, in the loop
    total_copies = initial_families + len(family_names) - len(named_starts)
    total_chromosomes = n_initial_chrom
    # each live driver paired with what it reads, resolved once: the names and the family ids are fixed
    live_rate_reads = [(src, live_target(src, named, module_map or {})) for src in live_rate_keys]
    live_ext_reads = [(src, live_target(src, named, module_map or {})) for src in live_ext_keys]

    # Recipient rules that read live gene content get their reader now, when the family ids and the
    # genomes exist. A family's own rule is keyed by the id the family was given, which is what each
    # of its genes carries.
    def with_reader(choice):
        rule, groups_c, traj_c = choice
        if isinstance(rule, Driven) and names_a_live_level(rule.driver):
            traj_c = _LiveOrderedContent(live_target(rule.driver, named, module_map or {}), gen, pos,
                                         counts)
        return rule, groups_c, traj_c

    transfer_to, group_of, to_traj = with_reader((transfer_to, group_of, to_traj))
    fam_choice = {named[declared[i].name]: with_reader(c) for i, c in fam_choice_prepared.items()}

    # eleven bare numbers on their default scopes — no modifier on any rate, none on any extent,
    # no per-lineage budget — is the common run, and it needs none of the loop's context machinery:
    # each total is scope(base) exactly, resolved here once. A rate whose base is None carries a
    # set_by, which is a modifier, so `plain` is False and its 0.0 is never read.
    plain = (not any_written and not any(r.modifiers for r in _rates.values()) and not any_per_lineage
             and not any(e.has_modifiers for e in (dup_ext, los_ext, tra_ext, inv_ext, trp_ext,
                                                   trl_ext)))
    dup_base, los_base, tra_base = dup.base or 0.0, los.base or 0.0, tra.base or 0.0
    org_base, inv_base, trp_base = org.base or 0.0, inv.base or 0.0, trp.base or 0.0
    trl_base, fis_base, fus_base = trl.base or 0.0, fis.base or 0.0, fus.base or 0.0
    cor_base, clo_base = cor.base or 0.0, clo.base or 0.0
    no_weights: dict = {}    # what `w` is when nothing is driven: read by .get, never written
    if plain:
        def _ext_ctx(k):
            # nothing in this run carries a modifier (`plain` pinned every extent to
            # has_modifiers=False), so an extent's sample() reads none of this: an empty context
            # is the same multiplication by 1.0, without rebuilding the loop's context per event.
            return {}

    # --- building a lineage's row ------------------------------------------------------------------
    # For bringing a row up to date from what changed in it: which live drivers read each family (its
    # presence, its module's completion), which read the gene count, and, per event class, which
    # declared families' own rates read each driver. All fixed for the run, so worked out once, in
    # lists, so the order a row visits them in is the run's and not a set's.
    drivers_of_family: dict[int, list] = collections.defaultdict(list)
    count_drivers: list = []
    for src, target in live_rate_reads:
        if target is None:
            count_drivers.append((src, target))
        else:
            for fam in (target if isinstance(target, tuple) else (target,)):
                drivers_of_family[fam].append((src, target))
    readers_of: dict[str, dict] = {key: collections.defaultdict(list) for key in own_keys}
    if any_written:
        for key in own_keys:
            for fam, gene_rate in fam_driven_by_id[key].items():
                for driver in gene_rate.reads:
                    readers_of[key][driver].append(fam)
    counts.recording = any_rows       # a row is brought up to date from the families that changed

    def run_unit_on(key, k, dk) -> float:
        """The run's rate per gene as lineage ``k`` reads it, for the families that did not write one
        of their own. The branch's multiplier belongs here and nowhere else in the own-rate
        arithmetic: it was written on the run's rate, so it scales the genes carrying that rate and
        not the genes whose family replaced it. Read through one function, so `own_sums` and the
        per-gene table `_own_table` hands an event cannot disagree about it."""
        unit = run_gene_rates[key].value(t, dk)
        return unit * lin_mult[key][alive[k]] if varying[key] else unit

    def w_from_rate_row(label, rate, k, size, n_chrom, values) -> float:
        """Lineage ``k``'s weight for one event class, read off the rate itself.

        Three shapes rather than one context built per call: a driven rate reads the branch's driver
        value, a rate varying among lineages reads the branch's drawn multiplier, and a rate doing
        both reads both. The driven-only shape is the expression it always was, so a conditioned run
        is unchanged to the last bit."""
        if not varying[label]:
            return rate.effective(copies=size, lineages=1, chromosomes=n_chrom, time=t,
                                  drivers=values)
        factor = lin_mult[label][alive[k]]
        if not driven[label]:
            return rate.effective(copies=size, lineages=1, chromosomes=n_chrom, time=t,
                                  carried_factor=factor)
        return rate.effective(copies=size, lineages=1, chromosomes=n_chrom, time=t,
                              drivers=values, carried_factor=factor)

    def w_from_family_row(label, k) -> float:
        """The same weight for a rate carrying a per-family draw **and** a per-lineage one: the
        family multipliers summed over the lineage's genes, times the branch's multiplier, times the
        run's unit rate. `fw` is left unscaled, because the segment pick reads it and a family that
        writes its own rate is measured against it (`set_own`)."""
        unit = _rates[label].effective(copies=1, lineages=1, chromosomes=1, time=t)
        return unit * lin_mult[label][alive[k]] * rows.fw[label][k]

    def family_total(label, rate, one, live=True) -> float:
        """One event class's total when the run draws per family: the unit rate times the weights
        summed over the live genes. A rate that also varies among lineages keeps that product per
        lineage, in its `w` row, because the branch's multiplier cannot be factored out of a sum over
        branches. Defined here rather than in the loop, which runs once per event."""
        if not live:
            return 0.0
        if label in rows.w:
            return rows.w[label].total
        return rate.effective(**one) * rows.fw[label].total

    def set_own(k: int, key: str, unit: float, owned: dict, written: float, covered: float,
                whole: float, largest: float) -> None:
        """Write lineage ``k``'s own-rate row for one event class. Every gene carries its family's own
        rate or the run's: the families that write one take ``covered`` of the lineage's ``whole``
        weight and give ``written``, and every other gene carries the run's rate on the rest."""
        if not owned:
            written = covered = 0.0       # nothing is left to carry a rounding error forward
        rows.own_sums[key].set(k, max(0.0, unit * (whole - covered) + written))
        rows.own_units[key][k] = unit
        rows.own_owned[key][k] = owned
        rows.own_written[key][k] = written
        rows.own_covered[key][k] = covered
        rows.own_largest[key][k] = largest

    def build_whole(k: int) -> None:
        """Build lineage ``k``'s row from its genome and counts alone."""
        counts.take_changed(k)            # a row built whole already holds whatever changed
        genome = gen[k]
        if any_driven:
            # A driven rate differs from lineage to lineage, so it is summed **over the living
            # lineages**, each read with its own driver value, its own gene count and its own
            # chromosome count — and the weights are kept, because the affected lineage must then be
            # drawn with them too. The gene count sits inside the weight, which is what makes a driven
            # per-copy rate a two-stage pick (a lineage, then a gene in it) rather than the one-stage
            # lineage draw a per-lineage rate takes. A per-family draw and a Driven cannot both be
            # set, so a `w` row built from a driver never sits beside an `fw` one; a `w` row built
            # from a per-lineage draw can, and `w_from_family_row` is where it is.
            values = {**{key: trajs[key].value(alive[k], t) for key in trajs},
                      **{src: _live_value(target, genome, counts.of(k))
                         for src, target in live_rate_reads}}
            rows.drivers[k] = values
        if any_driven or any_lineage:
            # the driver values this lineage reads, or None for a run with no driver at all: a rate
            # varying among lineages needs a row whether or not anything here is driven
            reads = rows.drivers[k] if any_driven else None
            size, n_chrom = _genome_size(genome), len(genome)
            for label in w_from_rate:
                rows.w[label].set(k, w_from_rate_row(label, _rates[label], k, size, n_chrom, reads))
        if any_family:
            # 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).
            for key, mult in fam_mult.items():
                rows.fw[key].set(k, sum(mult[g.family] for chrom in genome for g in chrom.genes))
            for label in w_from_family:
                rows.w[label].set(k, w_from_family_row(label, k))
        if any_written:
            # A family's own rate (SPEC §6): every gene carries a rate, its family's own when the
            # family writes one and the run's otherwise, and an event's total is their sum. The
            # acting lineage is drawn by its genes' summed rates and the segment by the mean rate of
            # the genes it covers, through the path the per-family draws take. The sum runs over the
            # families that write their own rate; every other gene carries the run's rate, so their
            # share is the lineage's whole weight less what the writing families take.
            held = counts.of(k)
            dk: dict[str, Any] = {"drivers": rows.drivers[k]} if any_driven else {}
            for key in own_keys:
                own_mult = fam_mult[key] if any_family else None
                unit = run_unit_on(key, k, dk)
                owned = {fam: value for fam, value in fam_fixed_by_id[key].items() if held[fam]}
                for fam, own_rate in fam_driven_by_id[key].items():
                    if held[fam]:
                        owned[fam] = own_rate.value(t, dk)
                written = covered = 0.0
                for fam, value in owned.items():
                    written += held[fam] * value
                    covered += held[fam] * (own_mult[fam] if own_mult is not None else 1.0)
                # the lineage's whole weight for this event: its summed per-family draws where the
                # run has them, its gene count otherwise
                whole = float(rows.fw[key][k] if any_family else _genome_size(genome))
                # every gene carries its family's own rate or the run's rate times its family's
                # draw, so the larger of those bounds all of them
                largest = max(unit * (own_mult.largest if own_mult is not None else 1.0),
                              max(owned.values(), default=0.0))
                set_own(k, key, unit, owned, written, covered, whole, largest)

    def bring_up_to_date(k: int) -> None:
        """Bring lineage ``k``'s row, built before, up to date from what changed in it since: the
        families whose copy number moved, the drivers those families move, and the families whose own
        rate reads a driver that moved. Nothing else in the row can have changed."""
        changed = counts.take_changed(k)
        genome = gen[k]
        held = counts.of(k)
        moved: list = []                  # the driver names whose value moved, in a fixed order
        if any_driven:
            values = rows.drivers[k]
            # the drivers that read a family that changed, and the gene count whenever the genome did
            asked = dict.fromkeys([pair for fam in changed for pair in drivers_of_family.get(fam, ())]
                                  + count_drivers)
            for src, target in asked:
                now = _live_value(target, genome, held)
                if now != values[src]:
                    if not moved:
                        values = dict(values)     # a new entry, so a kept snapshot is not edited
                    values[src] = now
                    moved.append(src)
            rows.drivers[k] = values
        if any_driven or any_lineage:
            # the driver values this lineage reads, or None for a run with no driver at all: a rate
            # varying among lineages needs a row whether or not anything here is driven
            reads = rows.drivers[k] if any_driven else None
            size, n_chrom = _genome_size(genome), len(genome)
            for label in w_from_rate:
                rows.w[label].set(k, w_from_rate_row(label, _rates[label], k, size, n_chrom, reads))
        if any_family and changed:
            for key, mult in fam_mult.items():
                weights = rows.fw[key]
                weights.set(k, weights[k] + sum(delta * mult[fam] for fam, delta in changed.items()))
            for label in w_from_family:
                rows.w[label].set(k, w_from_family_row(label, k))
        if any_written:
            dk: dict[str, Any] = {"drivers": rows.drivers[k]} if any_driven else {}
            size_now = float(_genome_size(genome))
            for key in own_keys:
                own_mult = fam_mult[key] if any_family else None
                fixed, driven_own, readers = fam_fixed_by_id[key], fam_driven_by_id[key], readers_of[key]
                unit = run_unit_on(key, k, dk)
                owned = rows.own_owned[key][k]
                written, covered = rows.own_written[key][k], rows.own_covered[key][k]
                largest = max(rows.own_largest[key][k],
                              unit * (own_mult.largest if own_mult is not None else 1.0))
                affected = [fam for fam in changed if fam in fixed or fam in driven_own]
                for src in moved:
                    affected.extend(readers.get(src, ()))
                if affected:
                    owned = dict(owned)
                    for fam in dict.fromkeys(affected):      # each family once, in a fixed order
                        weight = own_mult[fam] if own_mult is not None else 1.0
                        if fam in owned:                     # take out what it carried before
                            before = held[fam] - changed.get(fam, 0)
                            written -= before * owned.pop(fam)
                            covered -= before * weight
                        if held[fam]:                        # and put in what it carries now
                            value = fixed[fam] if fam in fixed else driven_own[fam].value(t, dk)
                            owned[fam] = value
                            written += held[fam] * value
                            covered += held[fam] * weight
                            largest = max(largest, value)
                whole = rows.fw[key][k] if any_family else size_now
                set_own(k, key, unit, owned, written, covered, whole, largest)

    bar = progress_bar(len(schedule), "genomes", unit="branch", enabled=progress)
    si = 0
    while si < len(schedule):
        if to_disk is not None:
            to_disk.flush_if_full()      # between steps, so a batch always holds whole events
        if _CHECK_COUNTS:
            _check_counts(gen, counts, rows)
        bar.to(si)
        n = total_copies
        k_alive = len(alive)
        c = total_chromosomes
        can_xfer = n > 0 and (k_alive >= 2 or self_transfer)
        next_species = schedule[si][0]
        # a family placed by `origins=` originates at a fixed instant, so it joins the horizon like
        # any other breakpoint: the waiting time can never step over it
        next_start = starts[start_i][0] if start_i < len(starts) else math.inf
        if plain:
            # no modifier on any rate or extent: each total is scope(base) exactly — a gene rate
            # times the live genes, a chromosome rate times the standing chromosomes, the two
            # originations per living lineage — none of them ever changes on its own (next_change
            # is inf), and nothing below reads a context or a weight.
            w = no_weights
            fw = None
            own_pick = {}
            r_dup = dup_base * n
            r_los = los_base * n
            r_tra = tra_base * n if can_xfer else 0.0
            r_inv = inv_base * n
            r_trp = trp_base * n
            r_trl = trl_base * n
            r_org = org_base * k_alive
            r_fis = fis_base * c
            r_fus = fus_base * c
            r_cor = cor_base * k_alive
            r_clo = clo_base * c
            horizon = min(next_species, next_start)
        else:
            ctx = {"copies": n, "lineages": k_alive, "chromosomes": total_chromosomes, "time": t}
            # A gene-level event counted PER LINEAGE is counted per lineage that HOLDS a gene: an
            # empty genome offers nothing to act on, so it must not take a share of the total and
            # then be picked with no victim inside it. Built only when some rate needs it, so the
            # per-copy path does exactly the work it did before.
            if any_per_lineage:
                gene_hosts = [k for k in range(k_alive) if _genome_size(gen[k])]
                gene_ctx = {**ctx, "lineages": len(gene_hosts)}
            else:
                gene_hosts, gene_ctx = None, ctx
            # Every number below reads one lineage on its own, so the step builds the rows of the
            # lineages an event changed and reads the rest as they stand (`_LineageRows`).
            if any_rows:
                for k, whole in rows.take_stale():
                    if whole:
                        build_whole(k)
                    else:
                        bring_up_to_date(k)
                if _CHECK_ROWS:
                    kept = rows.snapshot()
                    rows.touched_all()
                    for k, _whole in rows.take_stale():
                        build_whole(k)
                    rows.check_against(kept)
            w = rows.w
            fw = rows.fw if any_family else None

            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 w[label].total if label in w else pooled

            if fw is not None:   # the run draws per family: the same test as `any_family`
                one = {"copies": 1, "lineages": 1, "chromosomes": 1, "time": t}
                r_dup = family_total("duplication", dup, one, live=bool(n))
                r_los = family_total("loss", los, one, live=bool(n))
                r_tra = family_total("transfer", tra, one, live=can_xfer)
                r_inv = family_total("inversion", inv, one, live=bool(n))
                r_trp = family_total("transposition", trp, one, live=bool(n))
                r_trl = family_total("translocation", trl, one, live=bool(n))
            else:
                # each gene-level rate is read in the context its own scope asks for: `gene_ctx` counts
                # only the occupied genomes, which is what a per-lineage budget is counted over
                def _gc(label):
                    return gene_ctx if per_lineage[label] else ctx

                r_dup = _r("duplication", dup.effective(**_gc("duplication")) if n else 0.0, live=bool(n))
                r_los = _r("loss", los.effective(**_gc("loss")) if n else 0.0, live=bool(n))
                r_tra = _r("transfer", tra.effective(**_gc("transfer")) if can_xfer else 0.0,
                           live=can_xfer)
                r_inv = _r("inversion", inv.effective(**_gc("inversion")) if n else 0.0, live=bool(n))
                r_trp = _r("transposition", trp.effective(**_gc("transposition")) if n else 0.0,
                           live=bool(n))
                r_trl = _r("translocation", trl.effective(**_gc("translocation")) if n else 0.0,
                           live=bool(n))
            if any_written:
                if _CHECK_OWN_SUMS:
                    held_now = [counts.of(k) for k in range(k_alive)]
                    for key in own_keys:
                        _check_own_sums(own_pick[key][0], held_now, own_pick[key][1])
                if "duplication" in own_pick:
                    r_dup = own_pick["duplication"][0].total
                if "loss" in own_pick:
                    r_los = own_pick["loss"][0].total
                if "transfer" in own_pick:
                    r_tra = own_pick["transfer"][0].total if can_xfer else 0.0
            r_org = _r("origination", org.effective(**ctx))                 # per lineage
            r_fis = _r("fission", fis.effective(**ctx) if c else 0.0, live=bool(c))  # per chromosome
            r_fus = _r("fusion", fus.effective(**ctx) if c else 0.0, live=bool(c))
            r_cor = _r("chromosome_origination", cor.effective(**ctx))      # per lineage (de-novo replicon)
            r_clo = _r("chromosome_loss", clo.effective(**ctx) if c else 0.0, live=bool(c))
            # the next instant a rate changes on its own, asked only of the rates that ever do — a
            # family's own schedule among them (`timed_rates`)
            horizon = min(next_species, next_start)
            for timed in timed_rates:
                horizon = min(horizon, timed.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(k_alive)), default=math.inf))
            def _ext_ctx(k):
                """The context an extent is sampled in, on the lineage the event landed on.

                It cannot be built before the lineage is drawn, because a driven extent is read on
                the **acting** lineage at the instant the event fires — which is also why an extent
                adds no Gillespie breakpoint and never enters the horizon above (SPEC §6). With no
                driven extent this is the same context the rates were read in.

                The rest of `ctx` — the gene, lineage and chromosome counts — goes with it, because
                `Modifier.implemented_for` promises this engine supplies them and a modifier of
                your own is admitted onto an extent by the same gate that admits it onto a rate.
                Handing an extent a thinner context meant one gate certifying two different
                contracts: a modifier written the documented way read zeros, and one with a
                required keyword died mid-run."""
                # `ctx` was snapshotted at the top of the loop, before `t` advanced to the firing
                # instant, so `time` has to be taken fresh: an extent's own breakpoints are kept out
                # of the horizon, so a schedule's breakpoint routinely falls inside a stretch, and
                # reading the stale `t` would size the event on the wrong side of it.
                if not any_ext_driven:
                    return {**ctx, "time": t}
                return {**ctx, "time": t,
                        "drivers": {**{key: resolved[key].value(alive[k], t) for key in resolved},
                                    **{src: _live_value(target, gen[k], counts.of(k))
                                       for src, target in live_ext_reads}}}
        total = (r_dup + r_los + r_org + r_tra + r_inv + r_trp + r_trl
                 + r_fis + r_fus + r_cor + r_clo)

        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, rows.genes, fw, fam_mult, "duplication", dup_ext,
                                             _ext_ctx, w.get("duplication"),
                                             gene_hosts if per_lineage["duplication"] else None,
                                             own=own_pick.get("duplication"))
                    if picked is not None:
                        k, ci, j, m = picked
                        copied = _run_families(gen[k][ci], j, m)
                        if not _run_over_cap(counts.of(k), copied, cap):
                            total_copies += _duplicate(gen[k][ci], j, m, tree.nodes[alive[k]], t,
                                                       events, event_positions, new_gene)
                            counts.added_all(k, copied)
                            rows.touched(k, gen[k])
                elif r < b_los:
                    picked = _pick_event_run(rng, gen, n, rows.genes, fw, fam_mult, "loss", los_ext,
                                             _ext_ctx, w.get("loss"),
                                             gene_hosts if per_lineage["loss"] else None,
                                             own=own_pick.get("loss"))
                    if picked is not None:
                        k, ci, j, m = picked
                        taken = _run_families(gen[k][ci], j, m)
                        if _lose_at(gen[k][ci], j, m, tree.nodes[alive[k]], t, events, event_positions):
                            total_copies -= m
                            counts.removed_all(k, taken)
                            rows.touched(k, gen[k])
                elif r < b_org:
                    # origination is per lineage: a uniform lineage, or one drawn by its own rate
                    # when that rate is driven (the same weights the total was summed with)
                    k = (w["origination"].pick(rng) if "origination" in w
                         else int(rng.integers(k_alive)))
                    counts.added(k, _originate(gen[k], tree.nodes[alive[k]], t, events, event_positions,
                                               new_gene, new_family, rng))
                    total_copies += 1
                    rows.touched(k, gen[k])
                elif r < b_tra:
                    picked = _pick_event_run(rng, gen, n, rows.genes, fw, fam_mult, "transfer", tra_ext,
                                             _ext_ctx, w.get("transfer"),
                                             gene_hosts if per_lineage["transfer"] else None,
                                             own=own_pick.get("transfer"))
                    if picked is not None:                # driven: the weighted lineage is the DONOR
                        kd, cdi, jd, m = picked
                        delta, kr = _do_transfer(rng, tree, alive, gen, counts, kd, cdi, jd, m, t,
                                                 events, event_positions, new_gene, transfer_to,
                                                 replacement, self_transfer, depth, cap,
                                                 to_traj, group_of, fam_choice)
                        total_copies += delta
                        if kr is not None:   # the recipient's gene content changed, the donor's did not
                            rows.touched(kr, gen[kr])
                elif r < b_inv:
                    picked = _pick_event_run(rng, gen, n, rows.genes, fw, fam_mult, "inversion", inv_ext,
                                             _ext_ctx, w.get("inversion"),
                                             gene_hosts if per_lineage["inversion"] else None)
                    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, rows.genes, fw, fam_mult, "transposition", trp_ext,
                                             _ext_ctx, w.get("transposition"),
                                             gene_hosts if per_lineage["transposition"] else None)
                    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, rows.genes, fw, fam_mult, "translocation", trl_ext,
                                             _ext_ctx, w.get("translocation"),
                                             gene_hosts if per_lineage["translocation"] else None)
                    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:
                    picked = _pick_chromosome(rng, gen, c, rows.chromosomes, w.get("fission"))
                    if picked is not None:
                        k, ci = picked
                        dc, dg = _fission(gen[k], ci, tree.nodes[alive[k]], t, chromosome_events,
                                          new_chromosome, rng)
                        total_chromosomes += dc
                        total_copies += dg
                        if dc or dg:
                            rows.touched(k, gen[k])
                elif r < b_fus:
                    picked = _pick_chromosome(rng, gen, c, rows.chromosomes, w.get("fusion"))
                    if picked is not None:
                        k, ci = picked
                        dc, dg = _fusion(gen[k], ci, tree.nodes[alive[k]], t, chromosome_events,
                                         new_chromosome, rng)
                        total_chromosomes += dc
                        total_copies += dg
                        if dc or dg:
                            rows.touched(k, gen[k])
                elif r < b_cor:
                    # chromosome origination is per lineage, uniform or driven, exactly as origination
                    k = (w["chromosome_origination"].pick(rng)
                         if "chromosome_origination" in w else int(rng.integers(k_alive)))
                    dc, dg = _chromosome_originate(gen[k], tree.nodes[alive[k]], t, chromosome_events,
                                                   new_chromosome)
                    total_chromosomes += dc
                    total_copies += dg
                    if dc or dg:
                        rows.touched(k, gen[k])
                else:
                    picked = _pick_chromosome(rng, gen, c, rows.chromosomes, w.get("chromosome_loss"))
                    if picked is not None:
                        k, ci = picked
                        taken = [g.family for g in gen[k][ci].genes]
                        dc, dg = _chromosome_lose(gen[k], ci, tree.nodes[alive[k]], t, events,
                                                  event_positions, chromosome_events)
                        if dc:
                            counts.removed_all(k, taken)
                        total_chromosomes += dc
                        total_copies += dg
                        if dc or dg:
                            rows.touched(k, gen[k])
                continue

        if horizon == next_species:  # advance to the tree's next event(s); process the whole tie-batch
            t = next_species
            if time_varying:             # time moved, and something here moves with it
                rows.touched_all()
            while si < len(schedule) and schedule[si][0] == t:
                i = schedule[si][1]
                g = gen[pos[i]]
                if to_disk is None:
                    genomes[i] = tuple(Chromosome(c.id, c.topology, tuple(c.genes)) for c in g)  # freeze
                else:
                    to_disk.branch_ended(i, g)     # written now, and not kept
                total_copies -= sum(len(c.genes) for c in g)
                total_chromosomes -= len(g)
                inherited = counts.retired(pos[i])  # what the daughters below inherit, if any
                rows.retired(pos[i])
                retire(alive, gen, pos, pos[i])
                node = tree.nodes[i]
                if node.children:  # a speciation: re-mint every chromosome and gene id
                    child_genomes: dict[int, list[Chromosome]] = {c: [] for c in node.children}
                    for pchrom in g:
                        dcids = []
                        per_daughter: list[list[GeneEdge]] = []
                        for c in node.children:
                            dcid = new_chromosome()
                            dcids.append(dcid)
                            dgenes, edges = [], []
                            for old in pchrom.genes:  # ZOMBI1: the gene ends and continues, fresh id
                                ng = new_gene(old.family, old.strand)
                                dgenes.append(ng)
                                edges.append(GeneEdge(t, "speciation", c, old.family, ng.id,
                                                   parent=old.id))
                            per_daughter.append(edges)
                            child_genomes[c].append(Chromosome(dcid, pchrom.topology, dgenes))
                        # the ids are minted daughter by daughter (which is what fixes them), but a
                        # gene's two edges are recorded together: one gene ending is one event, and
                        # the log writes it as one row naming both daughters
                        for gene_edges in zip(*per_daughter):
                            events.extend(gene_edges)
                        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)
                        counts.entered_like(inherited)   # a re-id of the parent: same families
                        rows.entered(cg)
                        total_copies += sum(len(ch.genes) for ch in cg)
                        total_chromosomes += len(cg)
                si += 1
        elif start_i < len(starts) and horizon == next_start:
            # a placed family arrives — the ordinary origination event, at a time and on a lineage
            # that were chosen rather than drawn. The lineage is live by construction (its time was
            # checked against that branch's own life), and a tie with the tree's schedule falls to
            # the branch above, so the daughters have entered by the time this runs.
            t = horizon
            if time_varying:
                rows.touched_all()
            while start_i < len(starts) and starts[start_i][0] == t:
                _, lineage, fam = starts[start_i]
                _originate(gen[pos[lineage]], tree.nodes[lineage], t, events, event_positions,
                           new_gene, new_family, rng, family=fam)
                counts.added(pos[lineage], fam)
                rows.touched(pos[lineage], gen[pos[lineage]])
                total_copies += 1
                start_i += 1
        else:
            t = horizon  # a skyline breakpoint: advance and re-evaluate the (now changed) rate
            rows.touched_all()

    bar.close()
    if _CHECK_COUNTS:
        _check_counts(gen, counts, rows)
    links = links_of({**_rates, **_extents}, transfer_to, declared, fam_driven_rates, fam_transfer_to,
                     module_map or {})
    # a declared family's own rate, fixed or driven, replaces the run's, so no draw reaches it there
    own = {key: set(fam_fixed_by_id[key]) | set(fam_driven_by_id[key]) for key in own_keys}
    multipliers = multipliers_of(fam_mult, ORDERED_TARGETS, own) if any_family else {}
    per_lineage_multipliers = (lineage_multipliers_of(lin_mult, tree.labels(),
                                                     ORDERED_LINEAGE_TARGETS)
                               if any_lineage else {})
    if to_disk is not None:
        return to_disk.close(seed=seed, links=links, initial_genome=initial_genome, named=len(named),
                             family_multipliers=multipliers,
                             lineage_multipliers=per_lineage_multipliers)
    return OrderedGenomesResult(tree, genomes, events, rearrangements, chromosome_events, seed,
                                named, module_map, event_positions, initial_genome, links,
                                multipliers, per_lineage_multipliers)

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, deletion=0.0, deletion_extent=5.0, insertion=0.0, insertion_extent=5.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, modules=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 chromosome events. 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.
• deletion (per lineage, mean deletion_extent) removes an arc as an indel: the same material goes, but no copy lineage ends, so it is recorded in deletions and not in the genealogy, and its breakpoints do not cut the root partition. The pair divides like this — loss changes what a lineage has (a copy dies, a gene can go whole), deletion changes how much of a surviving copy it carries. In practice that is a difference of scale: hundreds to thousands of base pairs against ones to tens.
• insertion (per lineage, mean insertion_extent) lays down a run of novel spacer at a legal position — the twin of deletion, and origination without the gene. Novel DNA descends from nothing, so it arrives on a fresh source under a fresh copy lineage and is recorded in events as a root; the one breakpoint it makes to open a gap for itself is indel-made, so the block it landed inside is not cut in two. The pair divides the same way: origination brings a new gene family into the run, insertion brings sequence.
• 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", "distance" / a Distance, Clades({...}, Between({...})) or Recipients().weighted_by(driver, mapping) — see below; 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 rates are 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.)

Driving out of the genome. A gene declared by a GFF is named, and a named gene answers is it in this lineage, right now? — result.presence("dnaA"), the driver the other two resolutions already hand out. modules={"flagellum": ["flgA", "flgB"]} groups declared genes, so result.completion("flagellum") gives the fraction of the group a lineage carries; a module changes nothing about how the genome evolves.

Conditioning (a trait drives who receives). transfer_to = Recipients().weighted_by(driver, mapping) weights the candidate recipients by another level, and the numbers are weights, not rate multipliers: they are normalised across the candidates, so they leave the total amount of transfer alone and only redistribute it (SPEC §5, a weight, not a rate). Weight 0 means "cannot receive"; when every candidate weighs 0 the transfer does not happen. A Between({...}) mapping reads the donor's value too, so transfer can be steered between guilds rather than merely into one, and Clades({...}, Between({...})) is the same steering by named clade, read off the tree instead of a driver. Since a transfer here is additive, steering changes only which lineage the arc lands on — never how much DNA moves, and never anything about the donor.

Source code in zombi2/genomes/nucleotide.py

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,
                                deletion=0.0, deletion_extent=5.0,
                                insertion=0.0, insertion_extent=5.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, modules=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 chromosome events.
    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.
    - ``deletion`` (**per lineage**, mean ``deletion_extent``) removes an arc as an **indel**: the
      same material goes, but no copy lineage ends, so it is recorded in ``deletions`` and not in the
      genealogy, and its breakpoints do not cut the root partition. The pair divides like this —
      ``loss`` changes what a lineage *has* (a copy dies, a gene can go whole), ``deletion`` changes
      how much of a surviving copy it *carries*. In practice that is a difference of scale: hundreds
      to thousands of base pairs against ones to tens.
    - ``insertion`` (**per lineage**, mean ``insertion_extent``) lays down a run of **novel spacer**
      at a legal position — the twin of ``deletion``, and ``origination`` without the gene. Novel DNA
      descends from nothing, so it arrives on a fresh source under a fresh copy lineage and is
      recorded in ``events`` as a root; the one breakpoint it makes to open a gap for itself is
      indel-made, so the block it landed inside is not cut in two. The pair divides the same way:
      ``origination`` brings a new **gene family** into the run, ``insertion`` brings **sequence**.
    - ``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"``, ``"distance"`` / a `Distance`,
      ``Clades({...}, Between({...}))`` or ``Recipients().weighted_by(driver, mapping)`` — see below;
      ``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 rates
    are 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.)

    **Driving out of the genome.** A gene declared by a GFF is named, and a named gene answers *is it
    in this lineage, right now?* — ``result.presence("dnaA")``, the driver the other two resolutions
    already hand out. ``modules={"flagellum": ["flgA", "flgB"]}`` groups declared genes, so
    ``result.completion("flagellum")`` gives the fraction of the group a lineage carries; a module
    changes nothing about how the genome evolves.

    **Conditioning (a trait drives who receives).** ``transfer_to =
    Recipients().weighted_by(driver, mapping)``
    weights the candidate recipients by another level, and the numbers are **weights**, not rate
    multipliers: they are normalised across the candidates, so they leave the total amount of transfer
    alone and only redistribute it (SPEC §5, a weight, not a rate). Weight 0 means "cannot receive"; when
    every candidate weighs 0 the transfer does not happen. A ``Between({...})`` mapping reads the
    **donor's** value too, so transfer can be steered between guilds rather than merely into one, and
    ``Clades({...}, Between({...}))`` is the same steering by named clade, read off the tree instead
    of a driver. Since a transfer here is additive, steering changes only which lineage the arc lands
    on — never how much DNA moves, and never anything about the donor."""
    tree = as_tree(tree, level="genomes")
    # 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
    # chromosome rates. 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),
               ("deletion", deletion, PerLineage), ("insertion", insertion, 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}")
        rate = as_rate(spec, default_scope=want)
        # `rate.scope` holds the scope CLASS, not an instance — a scope constructor returns the rate
        # itself — so this is an identity test rather than an isinstance one.
        assert rate.scope is not None            # as_rate fills the level's default where none was written
        if rate.scope is not want:
            raise ValueError(
                f"{label} has a {rate.scope.__name__} scope, but the nucleotide engine reads {label} "
                f"as {want.__name__} and cannot read it any other way. Write {want.__name__}(...), "
                f"or drop the scope and let the level fill in its own.")
        for m in rate.modifiers:
            if isinstance(m, Driven):
                check_not_a_kernel(m.mapping, label=label)
            if not is_implemented(m, IMPLEMENTED_MODIFIERS, "genomes.nucleotide"):
                raise ValueError(
                    f"{label} carries {describe(m)}, which the nucleotide genome engine does not "
                    f"support. It takes {', '.join(cell_name(w) for w in IMPLEMENTED_MODIFIERS)}.")
        _rates[label] = rate
    def _as_bp_extent(spec, label, any_base=False):
        """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 `Geometric` for a **segmental** event — that 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 extent reaches the mutator as a *mean*, so a shape
        it cannot express would be accepted and then silently sampled as a geometric anyway.

        An **indel** takes any shape (``any_base``), because its cut set is unrestricted — every
        position is legal for one — so there is nothing left to re-weight against and its size can be
        drawn outright. That is what makes ``Fixed(1)`` mean one nucleotide rather than a geometric of
        mean one, and what puts a power law, the shape indel lengths actually take, within reach. The modifiers are the ones this resolution
        supports 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 any_base and not isinstance(e.base, Geometric):
            raise ValueError(
                f"{label} has a {type(e.base).__name__} base, but the nucleotide engine takes 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 isinstance(m, Driven):
                check_not_a_kernel(m.mapping, label=label)
            if not is_implemented(m, IMPLEMENTED_MODIFIERS, "genomes.nucleotide"):
                raise ValueError(
                    f"{label} carries {describe(m)}, which the nucleotide genome engine does not "
                    f"support — an extent takes the same modifiers a rate does here "
                    f"({', '.join(cell_name(w) for w in IMPLEMENTED_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")
    deletion_extent = _as_bp_extent(deletion_extent, "deletion_extent", any_base=True)
    insertion_extent = _as_bp_extent(insertion_extent, "insertion_extent", any_base=True)
    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,
                "deletion_extent": deletion_extent, "insertion_extent": insertion_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}")
    # who receives, validated in the one place all three resolutions share (SPEC §5): the mapping's
    # numbers are weights over the candidate recipients, never a rate multiplier
    transfer_to = resolve_transfer_to(transfer_to)
    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]
    # The one check that is about the run rather than about one argument: an extent too small to
    # cover a gene leaves the gene-level counters at ~0 however high the rates. Here, where the
    # layout is known and nothing has been drawn yet, and on the DECLARED genes only — a de-novo gene
    # arrives later at `origination_extent` and would drag the shortest one down to nothing.
    _warn_if_extents_cannot_reach_a_gene(specs, layouts, _rates, _extents)
    # Conditioning: a rate written with scaled_by 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 driver 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 undriven run is untouched.
    driven = {label: [m for m in r.modifiers if isinstance(m, Driven)] for label, r in _rates.items()}
    ext_driven = {label: [m for m in e.modifiers if isinstance(m, Driven)]
                  for label, e in _extents.items()}
    by_key: dict[object, "Driven"] = {}
    for mods in (*driven.values(), *ext_driven.values()):
        for m in mods:
            by_key.setdefault(m.key, m)
    resolved = {}
    if by_key:
        resolved = {key: resolve_driver(m.driver, tree, step=m.step, path=m.path, level="genomes.nucleotide")
                    for key, m in by_key.items()}
        # a mapping whose states never occur leaves every lineage on the default factor, so the run
        # would secretly be the undriven model — refuse it here, naming the driver
        for mods in (*driven.values(), *ext_driven.values()):
            for m in mods:
                label = m.driver if isinstance(m.driver, str) else f"<{type(m.driver).__name__}>"
                check_mapping_fires(m.mapping, resolved[m.key].states(), driver_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)
    # The transfer_to choice is prepared **after** `trajs` is fixed, for the same reason: a driven
    # transfer_to is a weight, not a rate, so its trajectory must not join `trajs` and start adding
    # horizon breakpoints. `resolved` doubles as the driver cache, so a trait that drives both a rate
    # and who receives is loaded once and read from one trajectory.
    group_of, to_traj = prepare_transfer_to(tree, transfer_to, resolved, level="genomes.nucleotide")

    # by keyword rather than by position: eleven rates and seven extents in one call is exactly the
    # shape where inserting a twelfth silently shifts everything after it
    rates = _Rates(**{k: _rates[k] for k in _rates},
                   inversion_extent=inversion_extent, translocation_extent=translocation_extent,
                   transposition_extent=transposition_extent, loss_extent=loss_extent,
                   deletion_extent=deletion_extent, insertion_extent=insertion_extent,
                   duplication_extent=duplication_extent,
                   transfer_extent=transfer_extent, origination_extent=origination_extent,
                   inversion_probability=inversion_probability)
    depth = mean_root_to_tip(tree)                       # timescale for Distance weighting

    rng, seed = stream("genomes", seed)     # own stream, and a drawn seed if none was given
    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] = []
    deletions: list[Deletion] = []                      # the indel log: extent-changing, no genealogy
    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), gene_spans))
        # `initial`, as the copy lineage below is: a replicon the run *starts* with is not something
        # it did, so counting `origination` in either log gives the de-novo ones alone
        chromosome_events.append(ChromosomeEvent(root.birth_time, "initial", root.id, (), (cid,)))
        events.append(Origination(root.birth_time, root.id, cid, cp, source, 0, length, kind="initial"))

    # a module groups declared genes, so it can only be checked once the layout has named them — the
    # same validation the other two resolutions run against `family_names`, here against the GFF's
    module_map = resolve_modules(modules, gene_names)

    # 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)

    # thirteen bare numbers on their stated scopes — no modifier on any rate, none on any extent —
    # is the common run, and it needs none of the loop's context machinery: each total is
    # scope(base) exactly, resolved here once. A rate whose base is None carries a set_by, which is
    # a modifier, so `plain` is False and its 0.0 is never read.
    plain = (not any(r.modifiers for r in _rates.values())
             and not any(e.has_modifiers for e in _extents.values()))
    _b = {label: (rate.base or 0.0) for label, rate in _rates.items()}
    inv_b, trl_b, trp_b = _b["inversion"], _b["translocation"], _b["transposition"]
    los_b, del_b, ins_b = _b["loss"], _b["deletion"], _b["insertion"]
    dup_b, tra_b, org_b = _b["duplication"], _b["transfer"], _b["origination"]
    fis_b, fus_b = _b["fission"], _b["fusion"]
    cor_b, clo_b = _b["chromosome_origination"], _b["chromosome_loss"]
    no_weights: dict = {}    # what `w` is when nothing is driven: read by .get, never written

    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
        next_species = schedule[si][0]
        if plain:
            # no modifier on any rate or extent: each total is scope(base) exactly — a gene event
            # per living lineage, a chromosome event per standing chromosome — none of them ever
            # changes on its own (next_change is inf), and nothing below reads a context or a
            # weight.
            w = no_weights
            r_inv = inv_b * nlin
            r_trl = trl_b * nlin
            r_trp = trp_b * nlin
            r_los = los_b * nlin
            r_del = del_b * nlin
            r_ins = ins_b * nlin
            r_dup = dup_b * nlin
            r_tra = tra_b * nlin if can_xfer else 0.0
            r_org = org_b * nlin
            r_fis = fis_b * count
            r_fus = fus_b * count
            r_cor = cor_b * nlin
            r_clo = clo_b * count
            horizon = next_species
        else:
            # 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 driver is
            # byte-identical to before.
            w = {}
            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_del = _r("deletion", rates.deletion.effective(**ctx))
            r_ins = _r("insertion", rates.insertion.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))
            # 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.deletion.next_change(t), rates.insertion.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))
        total = (r_inv + r_trl + r_trp + r_los + r_del + r_ins + r_dup + r_tra + r_org + r_fis
                 + r_fus + r_cor + r_clo)

        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.has_modifiers:
                return e.base.mean()
            # the same `ctx` the rates were read in, not a thinner one: an extent's modifiers pass
            # the gate that admits a rate's, so they are promised the same context (see `_ext_ctx`
            # in the ordered engine, which had the same hole).
            # `time` fresh, not `ctx`'s: `ctx` predates `t = t_ev`, and an extent is read at the
            # instant the event fires (see `_ext_ctx` in the ordered engine, same hole).
            return e.mean(**{**ctx, "time": t},
                          drivers={key: resolved[key].value(alive[k], t) for key in resolved})

        def _ext_draw(label, k):
            """One **drawn** size in bp, for an indel. The segmental events are parameterised by the
            mean (`_ext`) because they sample a far end out of the legal cut set; an indel has no
            such restriction, so it draws its size outright and any shape works —
            `Extent.sample()` scales the draw rather than the distribution's parameter, so the base
            still means what it says. At least 1 bp: a zero-length indel is not an event."""
            e = _extents[label]
            if not e.has_modifiers:
                return max(1, int(e.base.sample(rng)))
            return max(1, int(e.sample(rng, **{**ctx, "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 a per-chromosome one."""
            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
                # at deletion=0 this collapses onto b_los, so the ladder is numerically what it was
                # and no run written before indels existed moves
                b_del = b_los + r_del
                b_ins = b_del + r_ins
                b_dup = b_ins + 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_del:
                    k = _pick("deletion")
                    total_length += _do_deletion(gen[k], alive[k], t, _ext_draw("deletion_extent", k),
                                                 rng, deletions)
                elif r < b_ins:
                    k = _pick("insertion")
                    total_length += _do_insertion(gen[k], alive[k], t, _ext_draw("insertion_extent", k),
                                                  rng, events, new_source, new_copy)
                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,
                                                 to_traj, group_of)
                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:              # 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, module_map,
                                   gene_strands=gene_strands, initial_genome=initial_genome,
                                   initial_sequence=initial_sequence, deletions=deletions)

Results

One result type per resolution, each carrying the true history behind the dataset: the gene trees, the event log, and the genomes themselves.

zombi2.genomes.FamilyGenomesResult dataclass

FamilyGenomesResult(complete_tree: Tree, node_genomes: dict[int, tuple[GeneCopy, ...]], edges: list[GeneEdge], seed: int | None, family_names: dict[str, int] = dict(), modules: dict[str, tuple[str, ...]] = dict(), initial_genome: tuple[GeneCopy, ...] = (), max_family_size: int | None = None, links: 'tuple[Link, ...]' = (), family_multipliers: 'dict[int, dict[str, float | None]]' = dict(), lineage_multipliers: 'dict[int, dict[str, float]]' = dict())

What simulate_genomes_family returns: the complete_tree it ran on, the final node_genomes at every node (extant and extinct), the events log (the compact source of truth), and the seed. The observed dataset is the extant tips, genomes. The phyletic profiles are derived from those tips on access, and write materialises the chosen outputs to disk.

genomes property

genomes: dict[str, tuple[GeneCopy, ...]]

The observed dataset — the genome at each extant tip, keyed by the tip name the tree writes: n5.

Keyed by name because the only thing anyone does with this is join it to the tree, or to another level grown on that tree, and both name their tips. genomes.tsv on disk is keyed by name too, so what you get in Python and what you get from a file are the same dataset — they used to share no keys at all, and nothing said so. This is the trait level's TraitsResult.values for gene content.

node_genomes is the run's own record: every node, extant and extinct and internal alike, keyed by node id. Use that one to join against complete_tree.nodes or the event log.

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 profiles.

events cached property

events: list[Event]

The genome events — one per row of genome_events.tsv, the same objects the writer formats.

edges is the finer record this is grouped from: one entry per gene-tree edge, so a duplication is two of them and a transfer likewise. That is the shape a gene tree is built out of, and it used to be what this attribute returned — which meant counting duplications in Python gave twice the file's number, and a filter on kind == "transfer" matched everything here and nothing there. One word, one meaning: an event is what the log has a row for.

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 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.node_genomes[node_id])

completion

completion(name: str)

A module's completion as a conditioning driver — ModuleCompletion, a number in [0, 1]: the fraction of the module's families a lineage carries.

Read it with a Curve, the way any continuous driver is read; a threshold goes there rather than here (lambda f: 8.0 if f > 0.8 else 1.0). A joint run reads the same number while it runs, as "genomes:module:<name>".

Source code in zombi2/genomes/family.py

def completion(self, name: str):
    """A module's completion as a **conditioning driver** — `ModuleCompletion`, a number in
    ``[0, 1]``: the fraction of the module's families a lineage carries.

    Read it with a `Curve`, the way any continuous driver is read; a threshold goes there rather
    than here (``lambda f: 8.0 if f > 0.8 else 1.0``). A joint run reads the same number while
    it runs, as ``"genomes:module:<name>"``."""
    from .presence import ModuleCompletion
    if name not in self.modules:
        raise KeyError(f"no module {name!r}; declared modules are {sorted(self.modules)}")
    return ModuleCompletion(self, name)

presence

presence(name: str)

The named family's presence as a conditioning driver — GenePresence.

has_family answers for one node; this answers for every lineage at every instant, which is what a driven rate needs::

switch=PerLineage(0.1).scaled_by(g.presence("tox"), {"present": 5.0, "absent": 1.0})

Source code in zombi2/genomes/family.py

def presence(self, name: str):
    """The named family's presence as a **conditioning driver** — `GenePresence`.

    ``has_family`` answers for one node; this answers for every lineage at every instant, which
    is what a driven rate needs::

        switch=PerLineage(0.1).scaled_by(g.presence("tox"), {"present": 5.0, "absent": 1.0})
    """
    from .presence import GenePresence
    if name not in self.family_names:
        raise KeyError(f"no named family {name!r}; declared families are "
                       f"{sorted(self.family_names)}")
    return GenePresence(self, name)

has_family

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

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

Source code in zombi2/genomes/family.py

def has_family(self, node_id: int, name: str) -> bool:
    """Whether the named family ``name`` (declared via ``families=``) is present — has ≥ 1 copy — in
    the genome at ``node_id``. The presence a joint ``scaled_by("genomes:<name>", …)`` reads as its driver."""
    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.node_genomes[node_id])

summary

summary() -> dict

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

Every count here is one per event, which is also what genome_events.tsv is now one row of. They are counted from the GeneEdge objects, and those are one per gene-tree edge: a duplication, a transfer and a speciation each end one gene and start two, so counting edges inflates them exactly 2×. A duplication's two edges share the parent gene they descend from, and a gene ends at exactly one event, so distinct parents are the events.

loss counts every gene that died, which under replacement is more than the log's loss rows: a copy displaced by an arriving transfer has no row of its own — it is the second parent of that transfer_replacing row, because its death and the transfer are one event. It is a loss of the gene tree all the same, and that is what this counts.

The family counts are the other thing nobody could reconcile: gene_trees/ holds a file pair per family that ever existed, while the run's summary line counts the ones that survived, so "96 gene families" sat next to 213 files with nothing to explain the gap. Both numbers are here, named.

origination counts only the families that arose during the run: the initial genome is logged as origination at the root's own start time, so a bare count of that kind is de-novo arrivals plus initial_families, which is a number nobody asked for. They are separate here.

And the cap, which was invisible. When max_family_size binds it discards duplications and arriving transfers, so realised rates fall below the declared ones — so this reports which families are sitting at it, because that is the signal a reader can act on.

empty_genomes is the other end of the same story. There is no floor at this resolution: loss is counted per gene copy and the last copy is a copy like any other, so a high loss rate can strip a lineage of every gene it has. That is a real outcome of the model, not a failure — but it is invisible in the outputs, because a genome with no genes writes no row in profiles.tsv and leaves no gene tree for a sequence to run down. This is the number that says it happened, before a reader wonders why the matrix is short. (The ordered and nucleotide resolutions do have a floor, but it is a statement about what a chromosome is — a loss never takes a chromosome below its last gene — not a bound on genome size.)

Source code in zombi2/genomes/family.py

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

    **Every count here is one per event**, which is also what ``genome_events.tsv`` is now one row
    of. They are counted from the `GeneEdge` objects, and those are one per gene-tree *edge*: a
    duplication, a transfer and a speciation each end one gene and start two, so counting edges
    inflates them exactly 2×. A duplication's two edges share the ``parent`` gene they descend
    from, and a gene ends at exactly one event, so distinct parents *are* the events.

    ``loss`` counts every gene that died, which under ``replacement`` is more than the log's
    ``loss`` rows: a copy displaced by an arriving transfer has no row of its own — it is the
    second parent of that ``transfer_replacing`` row, because its death and the transfer are one
    event. It is a loss of the gene tree all the same, and that is what this counts.

    The family counts are the other thing nobody could reconcile: ``gene_trees/`` holds a file pair
    per family that ever existed, while the run's summary line counts the ones that survived, so
    "96 gene families" sat next to 213 files with nothing to explain the gap. Both numbers are
    here, named.

    ``origination`` counts only the families that arose *during* the run: the initial genome is
    logged as origination at the root's own start time, so a bare count of that kind is de-novo
    arrivals plus ``initial_families``, which is a number nobody asked for. They are separate here.

    And the cap, which was invisible. When ``max_family_size`` binds it discards duplications and
    arriving transfers, so realised rates fall below the declared ones — so this reports which
    families are sitting at it, because that is the signal a reader can act on.

    ``empty_genomes`` is the other end of the same story. There is **no floor** at this
    resolution: loss is counted per gene copy and the last copy is a copy like any other, so a
    high loss rate can strip a lineage of every gene it has. That is a real outcome of the model,
    not a failure — but it is invisible in the outputs, because a genome with no genes writes no
    row in ``profiles.tsv`` and leaves no gene tree for a sequence to run down. This is the
    number that says it happened, before a reader wonders why the matrix is short. (The ordered
    and nucleotide resolutions do have a floor, but it is a statement about what a chromosome is
    — a loss never takes a chromosome below its last gene — not a bound on genome size.)"""
    t0 = self.complete_tree.nodes[self.complete_tree.root].birth_time
    counted = event_counts(self.edges, t0)    # shared with the ordered and nucleotide summaries

    extant = list(self.complete_tree.extant_leaves())
    born = {e.family for e in self.edges}
    surviving = {c.family for i in extant for c in self.node_genomes.get(i, ())}
    genes_per_genome = [len(self.node_genomes.get(i, ())) for i in extant]
    cells = [collections.Counter(c.family for c in self.node_genomes.get(i, ())) for i in extant]
    copies = [n for cell in cells for n in cell.values()]

    cap = self.max_family_size
    at_cap = sorted({fam for cell in cells for fam, n in cell.items()
                     if cap is not None and n >= cap})
    return {
        "level": "genomes",
        "seed": self.seed,
        "resolution": "family",
        # one number per EVENT — the same thing a row of genome_events.tsv is
        "events": counted,
        "families": {"born": len(born), "surviving": len(surviving),
                     "died_out": len(born) - len(surviving),
                     "named": len(self.family_names)},
        "extant_genomes": len(extant),
        "empty_genomes": sum(1 for i in extant if not self.node_genomes.get(i, ())),
        "genes_per_genome": _stats(genes_per_genome),
        "copies_per_family_per_genome": _stats(copies),
        # the cap made visible. `families_at_cap` is what to look at: a family sitting at the
        # ceiling had events discarded, so its realised rates are below the ones you declared.
        "family_size_cap": {
            "cap": cap,
            "families_at_cap": len(at_cap),
            "cells_at_cap": sum(1 for cell in cells for n in cell.values()
                                if cap is not None and n >= cap),
            "family_ids_at_cap": at_cap},
    }

write

write(directory, outputs=('events', 'profiles', 'genomes', 'initial_genome', 'gene_trees', 'species_tree', 'summary', 'links', 'family_multipliers', 'lineage_multipliers'), *, flat: bool = False, compact: bool = False, bundle: 'bool | None' = None) -> 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 under gene_trees/, each family's true genealogy. A family with no surviving copy writes no _extant file.

• "species_tree" → species_complete.nwk, the tree the run evolved along. Written because a directory of gene trees with no species tree is not a dataset anyone can use: every one of these outputs is indexed by that tree's node labels, and the truth a gene tree is compared against is the species tree it grew inside. A run written from Python used to leave it out entirely, so the quickstart handed back gene trees with nothing to compare them to and said nothing about it.

• "links" → links.tsv, the links the run read from its own gene content, one row per link (see zombi2.genomes.links); the header alone when there are none.
• "family_multipliers" → family_multipliers.tsv, each family's drawn rate multipliers, one row per family (see zombi2.genomes.multipliers); the header alone when no rate varies among families.
• "lineage_multipliers" → lineage_multipliers.tsv, each species branch's drawn rate multipliers, one row per branch; the header alone when no rate varies among lineages. The gene trees are two files per family, so they get a subdirectory rather than burying the tables above; flat=True writes everything into directory instead. compact=True writes them as one file, gene_trees.tsv, one row per family: a run of thousands of families is then a handful of files, which is what a cluster that limits file counts needs.

Source code in zombi2/genomes/family.py

def write(self, directory, outputs=("events", "profiles", "genomes", "initial_genome",
                                    "gene_trees", "species_tree", "summary", "links",
                                    "family_multipliers", "lineage_multipliers"), *,
          flat: bool = False, compact: bool = False,
          bundle: "bool | None" = None) -> 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`` under ``gene_trees/``,
      each family's true genealogy. A family with no surviving copy writes no ``_extant`` file.

    - ``"species_tree"`` → ``species_complete.nwk``, the tree the run evolved along. Written
      because a directory of gene trees with no species tree is not a dataset anyone can use:
      every one of these outputs is *indexed by* that tree's node labels, and the truth a gene
      tree is compared against is the species tree it grew inside. A run written from Python
      used to leave it out entirely, so the quickstart handed back gene trees with nothing to
      compare them to and said nothing about it.
    - ``"links"`` → ``links.tsv``, the links the run read from its own gene content, one row per
      link (see `zombi2.genomes.links`); the header alone when there are none.
    - ``"family_multipliers"`` → ``family_multipliers.tsv``, each family's drawn rate
      multipliers, one row per family (see `zombi2.genomes.multipliers`); the header alone when no
      rate varies among families.
    - ``"lineage_multipliers"`` → ``lineage_multipliers.tsv``, each species branch's drawn rate
      multipliers, one row per branch; the header alone when no rate varies among lineages.
    The gene trees are two files per family, so they get a subdirectory rather than burying the
    tables above; ``flat=True`` writes everything into ``directory`` instead. ``compact=True``
    writes them as one file, ``gene_trees.tsv``, one row per family: a run of thousands of
    families is then a handful of files, which is what a cluster that limits file counts needs.
    """
    compact = compact_or_bundle(compact, bundle, "write")   # its 0.52.0 name
    # An unknown token used to write nothing and exit clean — silent data loss you discover
    # three pipeline steps later, when the next tool has no input. The other levels have always
    # raised; these two did not.
    if unknown := [o for o in outputs if o not in self.OUTPUTS]:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(self.OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    # a run's directory describes that run: clear the per-unit directories this write is
    # about to fill, so nothing from a previous run survives inside them (see fresh_dirs)
    fresh_dirs(d, ("gene_trees",), flat)
    names = self.complete_tree.labels()      # e<id> for a lineage that died; n<id> for the rest
    if "events" in outputs:
        (d / "genome_events.tsv").write_text(events_tsv(self.edges, names), encoding="utf-8")
    if "profiles" in outputs:
        (d / "profiles.tsv").write_text(self.profiles.to_tsv(), encoding="utf-8")
    if "genomes" in outputs:
        (d / "genomes.tsv").write_text(self._genomes_tsv(), encoding="utf-8")
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv(), encoding="utf-8")
    if "gene_trees" in outputs:
        if compact:
            write_gene_tree_table(self.gene_trees, d / COMPACT_FILES["gene_trees"], names)
        else:
            write_gene_trees(self.gene_trees, grouped_dir(d, "gene_trees", flat), names)
    if "species_tree" in outputs:
        (d / "species_complete.nwk").write_text(self.complete_tree.to_newick() + "\n",
                                                encoding="utf-8")
    if "summary" in outputs:
        write_summary(d / "genome_summary.json", self.summary())
    if "links" in outputs:
        (d / "links.tsv").write_text(links_tsv(self.links), encoding="utf-8")
    if "family_multipliers" in outputs:
        (d / "family_multipliers.tsv").write_text(
            multipliers_tsv(self.family_multipliers, FAMILY_TARGETS), encoding="utf-8")
    if "lineage_multipliers" in outputs:
        (d / "lineage_multipliers.tsv").write_text(
            lineage_multipliers_tsv(self.lineage_multipliers, self.complete_tree.labels()),
            encoding="utf-8")

zombi2.genomes.OrderedGenomesResult dataclass

OrderedGenomesResult(complete_tree: Tree, node_genomes: dict[int, tuple[Chromosome, ...]], edges: list[GeneEdge], rearrangements: list[Inversion | Transposition | Translocation], chromosome_events: list[ChromosomeEvent], seed: int | None, family_names: dict[str, int] = dict(), modules: dict[str, tuple[str, ...]] = dict(), event_positions: list[EventPosition] = list(), initial_genome: tuple[Chromosome, ...] = (), links: 'tuple[Link, ...]' = (), family_multipliers: 'dict[int, dict[str, float | None]]' = dict(), lineage_multipliers: 'dict[int, dict[str, float]]' = dict())

What simulate_genomes_ordered() returns: the complete_tree it ran on, the final genomes at every node as tuples of Chromosome\ s, the shared gene-genealogy events log, the rearrangements (inversions, transpositions and translocations) 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.

genomes property

genomes: dict[str, tuple[Chromosome, ...]]

The observed dataset — the genome at each extant tip, keyed by the tip name the tree writes: n5.

Keyed by name because the only thing anyone does with this is join it to the tree, or to another level grown on that tree, and both name their tips. The written output is keyed by name too, so what you get in Python and what you get from a file are the same dataset. This is the trait level's TraitsResult.values for gene content.

node_genomes is the run's own record: every node, extant and extinct and internal alike, keyed by node id. Use that one to join against complete_tree.nodes or the event log.

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 profiles.

events cached property

events: list[Event]

The genome events — one per row of genome_events.tsv, the same objects the writer formats.

edges is the finer record this is grouped from: one entry per gene-tree edge, so a duplication is two of them and a transfer likewise. That is the shape a gene tree is built out of, and it used to be what this attribute returned — which meant counting duplications in Python gave twice the file's number, and a filter on kind == "transfer" matched everything here and nothing there. One word, one meaning: an event is what the log has a row for.

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 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.node_genomes[node_id] for g in chrom.genes)

completion

completion(name: str)

A module's completion as a conditioning driver — ModuleCompletion, a number in [0, 1]: the fraction of the module's families a lineage carries.

Read it with a Curve, the way any continuous driver is read; a threshold goes there rather than here (lambda f: 8.0 if f > 0.8 else 1.0).

Source code in zombi2/genomes/ordered.py

def completion(self, name: str):
    """A module's completion as a **conditioning driver** — `ModuleCompletion`, a number in
    ``[0, 1]``: the fraction of the module's families a lineage carries.

    Read it with a `Curve`, the way any continuous driver is read; a threshold goes there rather
    than here (``lambda f: 8.0 if f > 0.8 else 1.0``)."""
    from .presence import ModuleCompletion
    if name not in self.modules:
        raise KeyError(f"no module {name!r}; declared modules are {sorted(self.modules)}")
    return ModuleCompletion(self, name)

presence

presence(name: str)

The named family's presence as a conditioning driver — GenePresence.

has_family answers for one node; this answers for every lineage at every instant, which is what a driven rate needs::

switch=PerLineage(0.1).scaled_by(g.presence("tox"), {"present": 5.0, "absent": 1.0})

Source code in zombi2/genomes/ordered.py

def presence(self, name: str):
    """The named family's presence as a **conditioning driver** — `GenePresence`.

    ``has_family`` answers for one node; this answers for every lineage at every instant, which
    is what a driven rate needs::

        switch=PerLineage(0.1).scaled_by(g.presence("tox"), {"present": 5.0, "absent": 1.0})
    """
    from .presence import GenePresence
    if name not in self.family_names:
        raise KeyError(f"no named family {name!r}; declared families are "
                       f"{sorted(self.family_names)}")
    return GenePresence(self, name)

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.node_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.node_genomes[node_id] for pos, g in enumerate(chrom.genes)]

write

write(directory, outputs=('events', 'profiles', 'gene_order', 'initial_genome', 'gene_trees', 'chromosome_events', 'species_tree', 'summary', 'links', 'family_multipliers', 'lineage_multipliers'), *, flat: bool = False, compact: bool = False, bundle: 'bool | None' = None) -> None

Materialise chosen outputs to directory (created if needed):

• "events" → two tables, because a run does two different things to a genome. genome_events.tsv is the gene genealogy — one row per event, in the format every resolution writes — with where each event happened beside it. rearrangement_events.tsv is the ancestry-neutral rearrangements: an inversion, a transposition or a translocation begins and ends no gene lineage, so it has no parents and no children and nothing to say in those columns. The two used to be one table, which meant nine columns empty on every rearrangement row and six on every genealogy row. Together with gene_order they are 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 under gene_trees/, each family's true genealogy — unchanged from the family resolution, position being orthogonal to it.

The gene trees are two files per family, so they get a subdirectory rather than burying the tables above; flat=True writes everything into directory instead. compact=True writes them as one file, gene_trees.tsv, one row per family. - "links" → links.tsv, the links the run read from its own gene content, one row per link (see zombi2.genomes.links); the header alone when there are none. - "family_multipliers" → family_multipliers.tsv, each family's drawn rate multipliers, one row per family (see zombi2.genomes.multipliers); the header alone when no rate varies among families.

Source code in zombi2/genomes/ordered.py

def write(self, directory, outputs=("events", "profiles", "gene_order", "initial_genome",
                                    "gene_trees", "chromosome_events", "species_tree",
                                    "summary", "links", "family_multipliers",
                                    "lineage_multipliers"), *,
          flat: bool = False, compact: bool = False,
          bundle: "bool | None" = None) -> None:
    """Materialise chosen ``outputs`` to ``directory`` (created if needed):

    - ``"events"`` → **two** tables, because a run does two different things to a genome.
      ``genome_events.tsv`` is the gene genealogy — one row per event, in the format every
      resolution writes — with **where** each event happened beside it.
      ``rearrangement_events.tsv`` is the ancestry-neutral rearrangements: an inversion, a
      transposition or a translocation begins and ends no gene lineage, so it has no parents and
      no children and nothing to say in those columns. The two used to be one table, which meant
      nine columns empty on every rearrangement row and six on every genealogy row. Together with
      ``gene_order`` they are 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`` under ``gene_trees/``,
      each family's true genealogy — unchanged from the family resolution, position being
      orthogonal to it.

    The gene trees are two files per family, so they get a subdirectory rather than burying the
    tables above; ``flat=True`` writes everything into ``directory`` instead. ``compact=True``
    writes them as one file, ``gene_trees.tsv``, one row per family.
    - ``"links"`` → ``links.tsv``, the links the run read from its own gene content, one row per
      link (see `zombi2.genomes.links`); the header alone when there are none.
    - ``"family_multipliers"`` → ``family_multipliers.tsv``, each family's drawn rate
      multipliers, one row per family (see `zombi2.genomes.multipliers`); the header alone when no
      rate varies among families.
    """
    compact = compact_or_bundle(compact, bundle, "write")   # its 0.52.0 name
    # An unknown token used to write nothing and exit clean — silent data loss you discover
    # three pipeline steps later, when the next tool has no input. The other levels have always
    # raised; these two did not.
    if unknown := [o for o in outputs if o not in self.OUTPUTS]:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(self.OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    # a run's directory describes that run: clear the per-unit directories this write is
    # about to fill, so nothing from a previous run survives inside them (see fresh_dirs)
    fresh_dirs(d, ("gene_trees",), flat)
    names = self.complete_tree.labels()   # e<id> for a lineage that died; n<id> for the rest
    if "events" in outputs:
        (d / "genome_events.tsv").write_text(
            _events_tsv(self.edges, self.event_positions, names), encoding="utf-8")
        (d / "rearrangement_events.tsv").write_text(
            rearrangement_events_tsv(self.rearrangements, names), encoding="utf-8")
    if "profiles" in outputs:
        (d / "profiles.tsv").write_text(self.profiles.to_tsv(), encoding="utf-8")
    if "gene_order" in outputs:
        (d / "gene_order.tsv").write_text(self._gene_order_tsv(names), encoding="utf-8")
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv(), encoding="utf-8")
    if "chromosome_events" in outputs:
        (d / "chromosome_events.tsv").write_text(
            chromosome_events_tsv(self.chromosome_events, self.complete_tree, names),
            encoding="utf-8")
    if "gene_trees" in outputs:
        if compact:
            write_gene_tree_table(self.gene_trees, d / COMPACT_FILES["gene_trees"], names)
        else:
            write_gene_trees(self.gene_trees, grouped_dir(d, "gene_trees", flat), names)
    if "species_tree" in outputs:            # the tree everything here is indexed by: without
        (d / "species_complete.nwk").write_text(   # it a directory of gene trees is not a dataset
            self.complete_tree.to_newick() + "\n", encoding="utf-8")
    if "summary" in outputs:
        write_summary(d / "genome_summary.json", self.summary())
    if "links" in outputs:
        (d / "links.tsv").write_text(links_tsv(self.links), encoding="utf-8")
    if "family_multipliers" in outputs:
        (d / "family_multipliers.tsv").write_text(
            multipliers_tsv(self.family_multipliers, ORDERED_TARGETS), encoding="utf-8")
    if "lineage_multipliers" in outputs:
        (d / "lineage_multipliers.tsv").write_text(
            lineage_multipliers_tsv(self.lineage_multipliers, self.complete_tree.labels(),
                                    ORDERED_LINEAGE_TARGETS), encoding="utf-8")

summary

summary() -> dict

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

The corrected event counts are the reason this exists. genome_events.tsv's loss rows undercount real losses whenever replacement is on, because a copy displaced by an arriving transfer has no row of its own; the migration guide names that as the change most likely to hand a returning user a plausible wrong number, and points them here. This file used to be written only at the family resolution — so the advice was sound and the remedy was absent at the two resolutions where the gap is larger (64% at ordered, measured).

event_counts is shared with the other two resolutions, so the three cannot drift. The rest is what this resolution has and the family core does not: where the genes sit, and what moved them.

Source code in zombi2/genomes/ordered.py

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

    The **corrected event counts** are the reason this exists. ``genome_events.tsv``'s ``loss``
    rows undercount real losses whenever ``replacement`` is on, because a copy displaced by an
    arriving transfer has no row of its own; the migration guide names that as the change most
    likely to hand a returning user a plausible wrong number, and points them here. This file used
    to be written only at the family resolution — so the advice was sound and the remedy was
    absent at the two resolutions where the gap is *larger* (64% at ordered, measured).

    `event_counts` is shared with the other two resolutions, so the three cannot drift. The rest
    is what this resolution has and the family core does not: where the genes sit, and what moved
    them."""
    t0 = self.complete_tree.nodes[self.complete_tree.root].birth_time
    extant = list(self.complete_tree.extant_leaves())
    return _ordered_summary(
        seed=self.seed, events=event_counts(self.edges, t0),
        born=len({e.family for e in self.edges}),
        surviving=len({g.family for i in extant for c in self.node_genomes.get(i, ())
                       for g in c.genes}),
        named=len(self.family_names), extant=len(extant),
        empty=sum(1 for i in extant if not any(c.genes for c in self.node_genomes.get(i, ()))),
        genes_per_genome=[sum(len(c.genes) for c in self.node_genomes.get(i, ())) for i in extant],
        chromosomes_per_genome=[len(self.node_genomes.get(i, ())) for i in extant],
        rearrangements=collections.Counter(type(r).__name__.lower() for r in self.rearrangements),
        chromosome_events=collections.Counter(e.kind for e in self.chromosome_events))

zombi2.genomes.NucleotideGenomesResult dataclass

NucleotideGenomesResult(complete_tree: Tree, node_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(), modules: dict[str, tuple[str, ...]] = dict(), gene_strands: dict[int, int] = dict(), initial_genome: NucleotideGenome = (lambda: NucleotideGenome([]))(), initial_sequence: dict[int, str] = dict(), deletions: list[Deletion] = list())

What 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, transposition), the chromosome_events (the chromosome network), and the seed. mosaic / trace_back / ancestry read a node's genome.

genomes property

genomes: dict[str, NucleotideGenome]

The observed dataset — the genome at each extant tip, keyed by the tip name the tree writes: n5.

Keyed by name because the only thing anyone does with this is join it to the tree, or to another level grown on that tree, and both name their tips. The written output is keyed by name too, so what you get in Python and what you get from a file are the same dataset. This is the trait level's TraitsResult.values for gene content.

node_genomes is the run's own record: every node, extant and extinct and internal alike, keyed by node id. Use that one to join against complete_tree.nodes or the event log.

indels property

indels: int

How many indels this run fired — deletions plus insertions. Both leave a lineage carrying part of a root block rather than all of it, which is the thing assembly() cannot yet put back together, so they are counted together and asked about in one place.

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 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 (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 root_blocks.

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 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 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 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.

genealogy property

genealogy: list[GeneEdge]

The run's genealogy as GeneEdge — the same table the family and ordered resolutions write, and what lands in genome_events.tsv.

A nucleotide run's own record, events, is interval-shaped: a copy lineage covers an extent, an event covers a sub-extent, and a duplication there mints a child without ending the parent (a split is not a birth). That is the right model for sequence, and the wrong shape for a gene tree. This is the translation onto the root-block partition, where a copy either covers a block in full or does not touch it — so a duplication is a bifurcation, writes two rows sharing a parent, and the ids are gene ids: the ones the gene trees, the alignments and the homology tables use. It is what the recovery already builds to derive gene_trees; writing it costs nothing extra.

describe

describe(node_id: int) -> str

One node's genome written out block by block, for reading by eye.

mosaic is the same information as data — (source, start, end, strand) per block — and this is it as text, one line per block, each labelled with the gene it is or intergene:

n2, chromosome 2
  [   0, 600) +  600 bp  intergene
  [ 600,1000) +  400 bp  gene 1

A block is a stretch of DNA with one unbroken ancestry, written as the interval it came from on the initial sequence and the strand it is read on, so a gene keeps its coordinates wherever it turns up and whichever way it points. A gene is always one block, because nothing may cut one; a run of intergene lines is spacer that events have cut apart.

Genes declared with a name (from a GFF) are named; the rest are numbered by family.

Source code in zombi2/genomes/nucleotide.py

def describe(self, node_id: int) -> str:
    """One node's genome written out block by block, for reading by eye.

    `mosaic` is the same information as data — ``(source, start, end, strand)`` per block — and
    this is it as text, one line per block, each labelled with the gene it is or ``intergene``:

        n2, chromosome 2
          [   0, 600) +  600 bp  intergene
          [ 600,1000) +  400 bp  gene 1

    A block is a stretch of DNA with one unbroken ancestry, written as the interval it came from
    on the **initial** sequence and the strand it is read on, so a gene keeps its coordinates
    wherever it turns up and whichever way it points. A gene is always one block, because
    nothing may cut one; a run of intergene lines is spacer that events have cut apart.

    Genes declared with a name (from a GFF) are named; the rest are numbered by family.
    """
    by_span = {span: fam for fam, span in self.gene_spans.items()}
    named = {fam: name for name, fam in self.gene_names.items()}
    # `node_label` and not an f-string: a lineage that died is `e5`, and the one place that
    # spelling is decided is the tree's own helper
    node = self.complete_tree.nodes.get(node_id)
    name = node_label(node_id, node.fate if node is not None else None)
    out = []
    for chromosome, blocks in self.mosaic(node_id).items():
        out.append(f"{name}, chromosome {chromosome}")
        for source, start, end, strand in blocks:
            fam = by_span.get((source, start, end))
            what = "intergene" if fam is None else named.get(fam, f"gene {fam}")
            out.append(f"  [{start:4d},{end:4d}) {'+' if strand == 1 else '−'} "
                       f"{end - start:4d} bp  {what}")
    return "\n".join(out)

completion

completion(name: str)

A module's completion as a conditioning driver — ModuleCompletion, a number in [0, 1]: the fraction of the module's genes a lineage carries.

Read it with a Curve, the way any continuous driver is read; a threshold goes there rather than here (lambda f: 8.0 if f > 0.8 else 1.0).

Source code in zombi2/genomes/nucleotide.py

def completion(self, name: str):
    """A module's completion as a **conditioning driver** — `ModuleCompletion`, a number in
    ``[0, 1]``: the fraction of the module's genes a lineage carries.

    Read it with a `Curve`, the way any continuous driver is read; a threshold goes there rather
    than here (``lambda f: 8.0 if f > 0.8 else 1.0``)."""
    from .presence import ModuleCompletion
    if name not in self.modules:
        raise KeyError(f"no module {name!r}; declared modules are {sorted(self.modules)}")
    return ModuleCompletion(self, name)

presence

presence(name: str)

The named gene's presence as a conditioning driver — GenePresence, the same reader the other two resolutions hand out, read off the gene's own recovered tree::

switch=PerLineage(0.1).scaled_by(g.presence("dnaA"), {"present": 5.0, "absent": 1.0})

A gene is named here by the GFF that declared it (its ID / Name); the evenly-spaced genes= layout lays its genes down unnamed.

Source code in zombi2/genomes/nucleotide.py

def presence(self, name: str):
    """The named gene's presence as a **conditioning driver** — `GenePresence`, the same reader
    the other two resolutions hand out, read off the gene's own recovered tree::

        switch=PerLineage(0.1).scaled_by(g.presence("dnaA"), {"present": 5.0, "absent": 1.0})

    A gene is named here by the GFF that declared it (its ``ID`` / ``Name``); the evenly-spaced
    ``genes=`` layout lays its genes down unnamed."""
    from .presence import GenePresence
    if name not in self.gene_names:
        raise KeyError(
            f"no named gene {name!r}; declared genes are {sorted(self.gene_names)}. A gene is "
            f"named by the GFF that declares it (its ID / Name attribute) — the evenly-spaced "
            f"genes= layout has no names to give.")
    return GenePresence(self, name)

block_of

block_of(family: int) -> int

The index in 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 `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, int, int]]]

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

The sub-range is what an indel makes necessary and is the whole of the presence bookkeeping: a deletion leaves a lineage holding a block minus a stretch of its middle, and an insertion opens a gap inside one, so a piece is a part of a block rather than all of it. Without indels every piece is the whole block — lo is 0 and hi its length — so the shape says the same thing it always did, at one extra pair of numbers.

To reconstruct a genome: take [lo, hi) of each piece's block 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. initial_assembly() does the same for the genome the run started with.

Every node votes on where the partition is cut (see _root_block_partition()), so a node's every event breakpoint is in it and no piece ever straddles two 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. Only an indel breakpoint, deliberately absent from the partition, makes a piece narrower than the block it indexes.

Source code in zombi2/genomes/nucleotide.py

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

    The sub-range is what an indel makes necessary and is the whole of the presence bookkeeping:
    a deletion leaves a lineage holding a block minus a stretch of its middle, and an insertion
    opens a gap inside one, so a piece is a *part* of a block rather than all of it. Without
    indels every piece is the whole block — ``lo`` is 0 and ``hi`` its length — so the shape says
    the same thing it always did, at one extra pair of numbers.

    To reconstruct a genome: take ``[lo, hi)`` of each piece's block 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.
    `initial_assembly()` does the same for the genome the run started with.

    Every node votes on where the partition is cut (see `_root_block_partition()`), so a node's
    every *event* breakpoint is in it and no piece ever straddles two 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. Only an indel
    breakpoint, deliberately absent from the partition, makes a piece narrower than the block it
    indexes."""
    tips = self._recover_blocks()[2]
    blocks = self.root_blocks
    what = node_label(node_id)
    out: dict[int, list[tuple[int, int, int, int, int]]] = {}
    for cid, pieces in self._pieces(self.node_genomes[node_id], what).items():
        named: list[tuple[int, int, int, int, int]] = []
        for (i, copy, strand, lo, hi) 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, lo, hi))
        out[cid] = named
    return out

initial_assembly

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

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

No gene id here, unlike 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 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, int, int]]]:
    """`assembly()` for `initial_genome`: ``{chromosome id: [(block, strand, lo, hi), …]}``.

    No gene id here, unlike `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 `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, lo, hi) for (i, _copy, strand, lo, hi) 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', 'species_tree', 'summary'), *, flat: bool = False, compact: bool = False, bundle: 'bool | None' = None) -> None

Materialise chosen outputs to directory (created if needed):

• "events" → three tables, because a nucleotide run records three different things. genome_events.tsv is the genealogy (genealogy) in the format every resolution writes, so one reader serves them all: one row per event, its participants named n<species>_g<copy>. block_events.tsv is this resolution's own record — the copy-lineage log over ancestral intervals, one row per interval an event touched, so an event spanning several blocks writes several rows sharing a time and kind. rearrangement_events.tsv is the ancestry-neutral rearrangements, which begin and end no lineage and so have nothing to put in parents and children. The last two are what read_nucleotide_genomes() replays.
• "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 under 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 under 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.

gene_trees, gff and bed are one file per family or per node — thousands, on a real genome times a real tree — so each gets a subdirectory rather than burying the tables above; flat=True writes everything into directory instead. compact=True writes the gene trees as one file, gene_trees.tsv, one row per family; gff and bed stay per node.

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", "species_tree", "summary"), *,
          flat: bool = False, compact: bool = False,
          bundle: "bool | None" = None) -> None:
    """Materialise chosen ``outputs`` to ``directory`` (created if needed):

    - ``"events"`` → **three** tables, because a nucleotide run records three different things.
      ``genome_events.tsv`` is the genealogy (`genealogy`) in the format *every* resolution
      writes, so one reader serves them all: one row per event, its participants named
      ``n<species>_g<copy>``. ``block_events.tsv`` is this resolution's own record — the
      copy-lineage log over **ancestral intervals**, one row per interval an event touched, so an
      event spanning several blocks writes several rows sharing a ``time`` and ``kind``.
      ``rearrangement_events.tsv`` is the ancestry-neutral rearrangements, which begin and end no
      lineage and so have nothing to put in ``parents`` and ``children``. The last two are what
      `read_nucleotide_genomes()` replays.
    - ``"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`` under ``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`` under ``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.

    ``gene_trees``, ``gff`` and ``bed`` are one file per family or per node — thousands, on a real
    genome times a real tree — so each gets a subdirectory rather than burying the tables above;
    ``flat=True`` writes everything into ``directory`` instead. ``compact=True`` writes the gene
    trees as one file, ``gene_trees.tsv``, one row per family; ``gff`` and ``bed`` stay per node.
    """
    compact = compact_or_bundle(compact, bundle, "write")   # its 0.52.0 name
    # An unknown token used to write nothing and exit clean — silent data loss you discover
    # three pipeline steps later, when the next tool has no input. The other levels have always
    # raised; these two did not.
    if unknown := [o for o in outputs if o not in self.OUTPUTS]:
        raise ValueError(f"unknown write outputs {unknown}; choose from {list(self.OUTPUTS)}")
    d = pathlib.Path(directory)
    d.mkdir(parents=True, exist_ok=True)
    # a run's directory describes that run: clear the per-unit directories this write is
    # about to fill, so nothing from a previous run survives inside them (see fresh_dirs)
    fresh_dirs(d, ("gene_trees", "gff", "bed"), flat)
    names = self.complete_tree.labels()   # e<id> for a lineage that died; n<id> for the rest
    if "events" in outputs:
        # Three tables, because they describe three things. `genome_events.tsv` is the genealogy
        # in the one format every resolution writes, so one reader serves them all;
        # `block_events.tsv` is this resolution's own interval record, which has no counterpart
        # elsewhere; `rearrangement_events.tsv` is what moved without ending anything. The first
        # two used to share a name, which made a nucleotide log look readable to the family reader
        # while meaning something else in the same columns; the third used to be rows inside the
        # second, with its copy and interval columns empty on every one of them.
        (d / "genome_events.tsv").write_text(events_tsv(self.genealogy, names), encoding="utf-8")
        (d / "block_events.tsv").write_text(
            _nucleotide_events_tsv([*self.events, *self.deletions], self.complete_tree, names),
            encoding="utf-8")
        (d / "rearrangement_events.tsv").write_text(
            rearrangement_events_tsv(self.rearrangements, names), encoding="utf-8")
    if "blocks" in outputs:
        (d / "blocks.tsv").write_text(self._blocks_tsv(names), encoding="utf-8")
    if "genes" in outputs:
        (d / "genes.tsv").write_text(self._genes_tsv(), encoding="utf-8")
    if "initial_genome" in outputs:
        (d / "initial_genome.tsv").write_text(self._initial_genome_tsv(), encoding="utf-8")
    if "chromosome_events" in outputs:
        (d / "chromosome_events.tsv").write_text(
            chromosome_events_tsv(self.chromosome_events, self.complete_tree, names),
            encoding="utf-8")
    if "gene_trees" in outputs:
        if compact:
            write_gene_tree_table(self.gene_trees, d / COMPACT_FILES["gene_trees"], names)
        else:
            write_gene_trees(self.gene_trees, grouped_dir(d, "gene_trees", flat), names)
    if "species_tree" in outputs:            # the tree everything here is indexed by: without
        (d / "species_complete.nwk").write_text(   # it a directory of gene trees is not a dataset
            self.complete_tree.to_newick() + "\n", encoding="utf-8")
    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:
            into = grouped_dir(d, token, flat)
            for label, genome in self._every_genome(names):
                (into / f"genome_{label}.{ext}").write_text(render(label, genome), encoding="utf-8")
    if "summary" in outputs:
        write_summary(d / "genome_summary.json", self.summary())

summary

summary() -> dict

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

The same correction as the other two resolutions: a GeneEdge is one gene-tree edge, so a duplication, a transfer and a speciation each end one gene and start two, and counting edges inflates them exactly 2×; this is the corrected count. (A transfer here is always additive, with no replacement option, so no displaced copy to fold in.) event_counts is shared with them, reading genealogy — the GeneEdge translation every resolution writes — so the three agree by construction rather than by three separate implementations agreeing by luck.

The unit here is the base pair, so there are no phyletic profiles and no per-family copy counts to report; what this resolution has instead is how much sequence there is and how it is divided.

Source code in zombi2/genomes/nucleotide.py

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

    The same correction as the other two resolutions: a `GeneEdge` is one gene-tree *edge*, so a
    duplication, a transfer and a speciation each end one gene and start two, and counting edges
    inflates them exactly 2×; this is the corrected count. (A transfer here is always additive,
    with no ``replacement`` option, so no displaced copy to fold in.) `event_counts` is shared
    with them, reading `genealogy` — the `GeneEdge` translation every resolution writes — so the
    three agree by construction rather than by three separate implementations agreeing by luck.

    The unit here is the base pair, so there are no phyletic profiles and no per-family copy
    counts to report; what this resolution has instead is how much sequence there is and how it is
    divided."""
    t0 = self.complete_tree.nodes[self.complete_tree.root].birth_time
    extant = list(self.complete_tree.extant_leaves())
    lengths = [self.node_genomes[i].length for i in extant if i in self.node_genomes]
    chrom_per_genome = [len(self.node_genomes[i].chromosomes) for i in extant if i in self.node_genomes]
    rearrangements = collections.Counter(type(r).__name__.lower() for r in self.rearrangements)
    chromosome = collections.Counter(e.kind for e in self.chromosome_events)
    # this resolution's own log, counted one number per EVENT — a row of `block_events.tsv` is
    # one ancestral *interval*, so an event spanning several blocks writes several rows. An
    # `Origination` splits the way `event_counts` splits it: the genome the run started with is
    # `initial`, and only a de-novo arrival is an `origination`.
    blocks = collections.Counter(
        (e.kind if isinstance(e, Origination) else type(e).__name__.lower())
        for e in self.events)
    return {
        "level": "genomes",
        "seed": self.seed,
        "resolution": "nucleotide",
        "events": event_counts(self.genealogy, t0),
        "extant_genomes": len(extant),
        # every gene the run holds, not only the ones declared at the start: `origination`
        # and `chromosome_origination` both mint genes into `gene_spans` as the run goes,
        # so a run that declared five can end with seven. It was called declared_genes.
        "genes": len(self.gene_spans),
        "base_pairs_per_genome": _stats(lengths),
        "chromosomes_per_genome": _stats(chrom_per_genome),
        # The same six kinds as `events`, so the two read side by side — and they must, because
        # an event here takes an arc of DNA while `events` counts gene-tree branchings, and
        # neither bounds the other: an arc covering three genes is three gene events, and an arc
        # covering none is zero. That is the whole answer to "why is duplication 0 when
        # block_events.tsv has thirteen rows" — see `_warn_if_extents_cannot_reach_a_gene`.
        "block_events": {k: blocks.get(k, 0)
                         for k in ("initial", "origination", "insertion", "duplication",
                                   "transfer", "loss", "speciation")},
        "rearrangements": {k: rearrangements.get(k, 0)
                           for k in ("inversion", "transposition", "translocation")},
        # the indel log, counted one number per EVENT like block_events above. It is neither a
        # genealogy event nor a rearrangement: no copy ends or begins, and material does go.
        "deletions": len(self.deletions),
        "base_pairs_deleted": sum(end - beg for e in self.deletions
                                  for (_cp, _src, beg, end) in e.deleted),
        "chromosome_events": dict(sorted(chromosome.items())),
    }

zombi2.genomes.StreamedRun dataclass

StreamedRun(directory: str, seed: 'int | None', n_families: int, n_events: int, outputs: tuple)

A genome run written straight to disk — what stream_to= returns, for a scale where a whole result would not fit in memory. At the family resolution the files are written family by family; at the ordered resolution they are written as the run goes. Thin by design: the outputs are the files and the disk is the handoff (the sequences level reads them back), so this carries where they are and how big the run was, not the run itself. n_events counts gene-tree edges, and n_families the families the run began.

path

path(output: str) -> str

The path of a written top-level file — e.g. path("events") → …/genome_events.tsv. Gene trees are not a single file; they live one pair per family under gene_trees/.

Source code in zombi2/genomes/_perfamily.py

def path(self, output: str) -> str:
    """The path of a written top-level file — e.g. ``path("events")`` → ``…/genome_events.tsv``.
    Gene trees are not a single file; they live one pair per family under ``gene_trees/``."""
    if output not in _STREAM_FILENAMES:
        raise KeyError(f"{output!r} is not a top-level streamed file (gene trees are under "
                       f"gene_trees/); files are {sorted(_STREAM_FILENAMES)}")
    return os.path.join(self.directory, _STREAM_FILENAMES[output])

Reading a run back

A written run is a genome run too: read_run reopens one from its directory, and the sequence level accepts a directory or a StreamedRun wherever it accepts a result.

zombi2.genomes.read_run

read_run(directory) -> FamilyGenomesResult

Reopen a genome run written to directory — the run object, from its files.

Takes a run directory in either layout (out/ with a genomes/ inside it, or a flat=True directory), or a StreamedRun handle. The event log is the source of truth, so gene_trees, events and the whole genealogy come back exactly; genomes and initial_genome come from their own tables when the run wrote them, and are empty otherwise.

This is what simulate_sequences calls when handed a path, so the two-command CLI pipeline and a two-step Python one are the same pipeline::

st = genomes.simulate_genomes_family(sp, ..., stream_to="run/")
sequences.simulate_sequences(st, model=hky85(), length=1000, seed=1)   # reads it back

The seed of a reopened run is the seed recorded on the handle when there is one, and None from a bare directory: the files are the run, and a number that did not produce them would be a lie. It is in the run's own run.zombi2 / genomes.log either way.

Source code in zombi2/genomes/read.py

def read_run(directory) -> FamilyGenomesResult:
    """Reopen a genome run written to ``directory`` — the run object, from its files.

    Takes a run directory in either layout (``out/`` with a ``genomes/`` inside it, or a ``flat=True``
    directory), or a `StreamedRun` handle. The event log is the source of truth, so ``gene_trees``,
    ``events`` and the whole genealogy come back exactly; ``genomes`` and ``initial_genome`` come from
    their own tables when the run wrote them, and are empty otherwise.

    This is what ``simulate_sequences`` calls when handed a path, so the two-command CLI pipeline and
    a two-step Python one are the same pipeline::

        st = genomes.simulate_genomes_family(sp, ..., stream_to="run/")
        sequences.simulate_sequences(st, model=hky85(), length=1000, seed=1)   # reads it back

    The ``seed`` of a reopened run is the seed recorded on the handle when there is one, and ``None``
    from a bare directory: the files are the run, and a number that did not produce them would be a
    lie. It is in the run's own ``run.zombi2`` / ``genomes.log`` either way.
    """
    seed = getattr(directory, "seed", None)                 # a StreamedRun carries its own
    directory = getattr(directory, "directory", directory)
    directory = os.fspath(directory)
    handoff = _resolve(directory)
    if os.path.exists(os.path.join(handoff, "blocks.tsv")):
        raise NotImplementedError(
            f"{handoff} is a nucleotide genome run (it has blocks.tsv), which reads back through "
            f"zombi2.genomes.nucleotide.read_nucleotide_genomes(directory, tree) — it needs the "
            f"blocks as well as the events, so it takes the tree explicitly.")
    tree = _tree(handoff, directory)
    with open(os.path.join(handoff, "genome_events.tsv"), encoding="utf-8") as f:
        edges = edges_from_tsv(f.read())
    initial: tuple[GeneCopy, ...] = ()
    initial_path = os.path.join(handoff, "initial_genome.tsv")
    if os.path.exists(initial_path):
        # read by column name: a family run writes `family copy`, and an ordered run writes the same two
        # among its layout columns (`chromosome topology position strand family copy`)
        with open(initial_path, encoding="utf-8") as f:
            header, *rows = [line.rstrip("\n").split("\t") for line in f if line.strip()]
        at_family, at_copy = header.index("family"), header.index("copy")
        initial = tuple(GeneCopy(gene_from_label(row[at_copy]), int(row[at_family])) for row in rows)
    links: tuple[Link, ...] = ()
    links_path = os.path.join(handoff, "links.tsv")
    if os.path.exists(links_path):
        with open(links_path, encoding="utf-8") as f:
            links = links_from_tsv(f.read())
    multipliers: dict[int, dict[str, "float | None"]] = {}
    multipliers_path = os.path.join(handoff, "family_multipliers.tsv")
    if os.path.exists(multipliers_path):
        with open(multipliers_path, encoding="utf-8") as f:
            multipliers = multipliers_from_tsv(f.read())
    # the other axis: a multiplier per species branch. Keyed back to node ids through the tree's own
    # labels, which the file writes and this directory ships, so the table lines up with the rest of
    # the result rather than carrying names nothing else here uses.
    per_lineage: dict[int, dict[str, float]] = {}
    lineage_path = os.path.join(handoff, "lineage_multipliers.tsv")
    if os.path.exists(lineage_path):
        with open(lineage_path, encoding="utf-8") as f:
            # node ids, because the labels are given — the `int | str` in the reader's own type is
            # the other call, which reads a table found without the tree it was written beside
            per_lineage = cast("dict[int, dict[str, float]]",
                               lineage_multipliers_from_tsv(f.read(), tree.labels()))
    return FamilyGenomesResult(complete_tree=tree, node_genomes=_genomes(handoff, tree), edges=edges,
                               seed=seed, initial_genome=initial, links=links,
                               family_multipliers=multipliers, lineage_multipliers=per_lineage)

Gene trees

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 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, names: 'dict[int, str] | None' = None, labels: str = 'events') -> str | None

Newick of the "extant" (default) or "complete" tree; None if it is empty. Leaves are n<species>_g<copy> — the copy and the branch it sits on, the same name the alignment FASTA records and the homology tables use, so a tip needs no translation to say which genome it came from. With annotate internal nodes carry <kind>_n<species>; branch lengths are time differences.

names is the run's node names (Tree.labels()), which is how a tree writes e<species> for a gene sitting on a lineage that died. The extant tree needs them too, which is not obvious: its leaves are all in living species, but a transfer node sits on the donor's branch, and a transfer out of a lineage that later died survives into the extant tree whenever the copy that moved has extant descendants. Without names that node was written n<species> for a species the extant tree does not contain and the complete tree calls e<species> — a label that resolved in neither file, on exactly the transfers whose donor is invisible.

labels picks what internal nodes are called. "events" (default) is the annotated form above. "copies" names every internal node the way the leaves are named, n<species>_g<copy> — the key the alignment and ancestral tables use — so a per-copy quantity (an ancestral sequence, a composition) can be attached to every branch of the tree, not only its tips. Event kinds are then not in the labels; they remain on the GeneNode objects.

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,
              names: "dict[int, str] | None" = None,
              labels: str = "events") -> str | None:
    """Newick of the ``"extant"`` (default) or ``"complete"`` tree; ``None`` if it is empty.
    Leaves are ``n<species>_g<copy>`` — the copy and the branch it sits on, the same name the
    alignment FASTA records and the homology tables use, so a tip needs no translation to say
    which genome it came from. With ``annotate`` internal nodes carry ``<kind>_n<species>``;
    branch lengths are time differences.

    ``names`` is the run's node names (`Tree.labels()`), which is how a tree writes
    ``e<species>`` for a gene sitting on a lineage that died. The **extant** tree needs them too,
    which is not obvious: its *leaves* are all in living species, but a transfer node sits on the
    **donor's** branch, and a transfer out of a lineage that later died survives into the extant
    tree whenever the copy that moved has extant descendants. Without ``names`` that node was
    written ``n<species>`` for a species the extant tree does not contain and the complete tree
    calls ``e<species>`` — a label that resolved in neither file, on exactly the transfers whose
    donor is invisible.

    ``labels`` picks what internal nodes are called. ``"events"`` (default) is the
    annotated form above. ``"copies"`` names every internal node the way the leaves are
    named, ``n<species>_g<copy>`` — the key the alignment and ancestral tables use — so a
    per-copy quantity (an ancestral sequence, a composition) can be attached to every
    branch of the tree, not only its tips. Event kinds are then not in the labels; they
    remain on the `GeneNode` objects.

    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."""
    if labels not in ("events", "copies"):
        raise ValueError(f"labels must be 'events' or 'copies', not {labels!r}")
    root = self.extant if which == "extant" else self.complete
    if root is None:
        return None
    return _to_newick(root, annotate, self.origination, names, labels) + ";"

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).

Who receives a transfer

transfer_to chooses who receives a horizontal transfer. It redistributes the transfers without changing how many happen. A weight of 0 means "cannot receive", and when every candidate weighs 0 the transfer does not fire. The numbers are weights, not rate multipliers, so a driven rule is written from its own entry point, transfer_to = Recipients().weighted_by(driver, mapping), with no base in front of it — a rate written there is an error.

zombi2.genomes.Distance dataclass

Distance(decay: float = 1.0)

A transfer_to weighting by relatedness: a recipient at patristic distance d from the donor gets weight exp(-decay × d / depth), where depth is the tree's mean root-to-tip time — so decay is scale-free (in units of tree depth), meaning the same across trees of different absolute timescales. transfer_to="distance" is Distance(decay=1.0).

zombi2.genomes.Clades dataclass

Clades(groups: dict, between: object)

A transfer_to weighting by named clades — the topological, donor-conditioned sibling of Distance. Each group is a clade of the species tree, and a Between kernel weights a candidate recipient by the pair (donor's clade, recipient's clade), so a transfer can be steered to run between two clades rather than within them — which the per-recipient weight of a Driven cannot express::

transfer_to = Clades({"A": ["n12", "n27"], "B": 40},
                     Between({("A", "B"): 1.0, ("B", "A"): 1.0}, default=0.0))

A clade is named either by a set of tips (a list — the clade is the subtree below their MRCA) or by a single node id (an int, or an "n<id>" label — the clade is that node's whole subtree). Groups must be disjoint; a lineage in none of them is in the implicit group "rest", usable as a kernel key. Membership is read from the tree (a clade is a fact about the tree, not another level), so this is a topological rule like "distance", resolved once per run — not a weighted_by reading another level, and needing no driver file.

zombi2.genomes.Between

Between(per_pair, default: float = 1.0)

A weight over ordered (donor-group, recipient-group) pairs — the 2-D kernel of the transfer choice of who receives (SPEC §5), the donor-conditioned sibling of Table::

Between({("A", "B"): 1.0, ("B", "A"): 1.0}, default=0.0)   # A↔B only, nothing else receives
Between({("A", "B"): 3.0})                                 # A→B 3× baseline, every other pair 1×

A Table weights a candidate recipient by that candidate's state alone; a Between weights it by the pair — the donor's group and the recipient's — which is what lets a transfer be steered to run between two groups rather than within them. It is therefore not a Mapping (a Mapping.multiplier reads one value): its weight() reads two, and the engine passes both. It is used in transfer_to — on its own as the kernel of a Clades rule (groups from the tree), or as the mapping of a Recipients().weighted_by(...) (groups from a trait). It is not a rate multiplier: a rate has no donor to condition on, so a Between on a rate is refused.

Keys are (from_group, to_group) pairs matched by string form, exactly like Table's states, so an integer-labelled group still finds its entry. default (1.0) is the weight for any pair not named — default=0.0 gives the "only the flows I name can happen" idiom, reusing the rule that a weight of 0 means the donor cannot send to that recipient group; when every candidate weighs 0 the transfer has nowhere to land and does not fire.

Source code in zombi2/params/mapping.py

def __init__(self, per_pair, default: float = 1.0) -> None:
    if not isinstance(per_pair, dict) or not per_pair:
        raise ValueError(
            f"Between needs a non-empty {{(from_group, to_group): weight}} dict, got {per_pair!r}")
    table: dict[tuple[str, str], float] = {}
    for pair, weight in per_pair.items():
        if not (isinstance(pair, tuple) and len(pair) == 2):
            raise ValueError(
                f"Between keys are (from_group, to_group) pairs, got {pair!r} — write "
                f"Between({{('A', 'B'): 1.0}}), the donor group first, the recipient group second")
        key = (str(pair[0]), str(pair[1]))  # groups matched by string form, like Table's states
        if key in table:
            raise ValueError(
                f"Between pairs collide as strings: {pair!r} and an earlier key both map to {key!r}")
        table[key] = _check_factor(weight, f"Between weight for {pair!r}")
    self.per_pair = table
    self.default = _check_factor(default, "Between default")

weight

weight(from_group: object, to_group: object) -> float

The weight for a transfer from a from_group donor to a to_group recipient — the named pair's weight, or default if the pair is unnamed.

Source code in zombi2/params/mapping.py

def weight(self, from_group: object, to_group: object) -> float:
    """The weight for a transfer from a ``from_group`` donor to a ``to_group`` recipient — the
    named pair's weight, or `default` if the pair is unnamed."""
    return self.per_pair.get((str(from_group), str(to_group)), self.default)

groups

groups() -> set

Every group named on either side of a pair — what a fires-check tests against the groups that actually occur, so a kernel naming only absent groups (a typo) can be caught.

Source code in zombi2/params/mapping.py

def groups(self) -> set:
    """Every group named on either side of a pair — what a fires-check tests against the groups
    that actually occur, so a kernel naming only absent groups (a typo) can be caught."""
    return {g for pair in self.per_pair for g in pair}