Instances and results#

A MaxMPDInstance holds the trees, tree weights, taxon costs and survival probabilities, parsed and validated once. Its methods optimize for one or many budgets, score taxon sets, and report the structure that determines the running time.

class forestpd.MaxMPDInstance(trees, costs=None, weights=None, survival=None, *, missing_branch_length='error', rooted=False)[source]#

A MaxMPD instance: phylogenetic trees on a shared taxon set, tree weights, taxon costs and survival probabilities.

The objective of a protected set S is the expected weighted PD of the taxa that survive: protected taxa survive surely, and every other taxon x independently with probability a(x). With all a(x) = 0 (the default) it is the plain weighted PD of S.

The trees are parsed and validated once; optimize() and evaluate_pd() can then be called any number of times.

Parameters:

  • trees (str | os.PathLike | Sequence) – The path of a file holding one Newick tree per line, or a list of Newick strings and/or PhyloZoo SemiDirectedPhyNetwork / DirectedPhyNetwork trees (with the phylozoo extra). Every Newick tree ends with ;. All trees must carry the same taxon set.

  • costs (Mapping[str, int] | None, optional) – Integer cost per taxon label, at least 1. None means unit costs; missing labels cost 1.

  • weights (Sequence[float] | Mapping[int, float] | None, optional) – Per-tree weight w_j, in the order of the trees (the list, or the file), or keyed by tree index. None means all 1.0.

  • survival (Mapping[str, float] | None, optional) – Probability a(x) in [0, 1) that taxon x survives when it is not protected. None means 0 for every taxon; missing labels have 0.

  • missing_branch_length (str | float, default "error") – Policy for edges without a branch length: "error", "zero", or a numeric default.

  • rooted (bool, default False) – Use rooted PD: the PD of a set includes the path from its taxa to the root of each tree, and the stem above the root (the root’s own branch length) when the set is non-empty. By default PD is unrooted and the rooting of the input is ignored.

Variables:

  • trees (tuple[PhyloTree, ...]) – The trees, with leaf id == taxon id. With rooted=True they carry one extra leaf at the root, which never appears in results.

  • registry (TaxonRegistry) – Taxon label <-> id mapping.

  • costs (numpy.ndarray) – int64[n]; cost per taxon id (by registry; with rooted=True this includes the hidden root leaf, so pair it with registry.labels, not with taxa).

  • weights (numpy.ndarray) – float64[t]; weight per tree.

  • survival (numpy.ndarray) – float64[n]; survival probability a(x) per taxon id.

Raises:

ForestPDInputError – If any input is malformed, or the trees do not share one taxon set.

Warns:

ForestPDWarning – If the trees look rooted (a [&R] comment, a PhyloZoo DirectedPhyNetwork, or a binary root in every tree) but rooted is False, or a tree is marked [&U] while rooted is True.

Examples

>>> import forestpd as fpd
>>> inst = fpd.MaxMPDInstance(["(A:1,B:1,(C:1,D:1):2);", "(A:1,C:1,(B:1,D:1):2);"])
>>> inst.optimize(budget=2).value
8.0
>>> inst.evaluate_pd({"A", "B"})
6.0
>>> rooted = fpd.MaxMPDInstance(["((A:1,B:1):1,(C:1,D:1):1);"], rooted=True)
>>> rooted.evaluate_pd({"A", "B"}), rooted.evaluate_pd({"A"})   # paths to the root count
(3.0, 2.0)
>>> risky = fpd.MaxMPDInstance(["(A:1,B:1,(C:1,D:1):2);"], survival={"C": 0.5, "D": 0.5})
>>> risky.evaluate_pd({"A", "B"})    # 1 + 1, plus 2 * 3/4 (C or D survives) + 1/2 + 1/2
4.5
classmethod from_phylotrees(trees, registry, costs=None, weights=None, survival=None)[source]#

Build an instance from trees already numbered against registry.

Parameters:

  • trees (Sequence[PhyloTree]) – The trees, with leaf id == taxon id of registry.

  • registry (TaxonRegistry) – The taxon registry the leaf ids refer to.

  • costs (Mapping[str, int] | Sequence[int] | None, optional) – Cost per taxon label, or a sequence of costs indexed by taxon id. None means unit costs.

  • weights (Sequence[float] | Mapping[int, float] | None, optional) – Per-tree weight w_j. None means all 1.0.

  • survival (Mapping[str, float] | Sequence[float] | None, optional) – Survival probability per taxon label, or a sequence indexed by taxon id. None means all 0.

Returns:

The instance.

Return type:

MaxMPDInstance

Raises:

ForestPDInputError – If a cost, weight or survival probability is invalid, the trees do not share one taxon set, or a tree has an internal vertex of degree 2 (suppress it first with suppress_degree_two()).

property n_trees#

Number of trees.

property rooted#

Whether PD is rooted.

property n_taxa#

Number of taxa.

property taxa#

Taxon labels, in id order.

optimize(budget, *, return_taxa=False, maf=None)[source]#

Find a taxon set S of cost at most budget maximizing the objective.

The objective is that of evaluate_pd().

Parameters:

  • budget (int) – Total integer budget B (>= 0).

  • return_taxa (bool, default False) – Also reconstruct an optimal taxon set.

  • maf (MAFComputer | None, optional) – Agreement-forest stage. None (default) uses HeuristicMAF.

Returns:

The optimal value, optionally an optimal taxon set, and diagnostics.

Return type:

MaxMPDResult

Raises:

  • ForestPDInputError – If budget is not a non-negative integer.

  • ForestPDInputError – If the MAF stage does not return an agreement forest (a partition of each cluster’s leaves), or maf is not a MAF computer, or a forest has more than 63 components.

  • ForestPDValueError – From a MAF stage that cannot handle a cluster (e.g. maffe and a polytomy).

optimize_many(budgets, *, return_taxa=False, maf=None)[source]#

optimize() for several budgets at once, from a single solve.

The DP tables computed for the largest budget contain the answer for every smaller one, so this costs one optimize() call plus one witness reconstruction per budget. When the budgets are large, the complementary tables (DP’, on the budget left unspent) are cheaper, or both tables split between small and large budgets (meta["tables"]).

Parameters:

  • budgets (Iterable[int]) – The budgets (each >= 0).

  • return_taxa (bool, default False) – Also reconstruct an optimal taxon set for each budget.

  • maf (MAFComputer | None, optional) – Agreement-forest stage. None (default) uses HeuristicMAF.

Returns:

One result per budget, in the given order.

Return type:

list[MaxMPDResult]

Raises:

  • ForestPDInputError – If a budget is not a non-negative integer.

  • ForestPDInputError – If the MAF stage does not return an agreement forest (a partition of each cluster’s leaves), or maf is not a MAF computer, or a forest has more than 63 components.

  • ForestPDValueError – From a MAF stage that cannot handle a cluster (e.g. maffe and a polytomy).

Examples

>>> import forestpd as fpd
>>> inst = fpd.MaxMPDInstance(["(A:1,B:1,(C:1,D:1):2);", "(A:1,C:1,(B:1,D:1):2);"])
>>> [r.value for r in inst.optimize_many([1, 2, 3])]
[0.0, 8.0, 10.0]
structure(maf=None)[source]#

The clusters and agreement-forest sizes, without solving.

Runs the first half of optimize(): the cluster reduction (cached on the instance, so a later optimize() does not repeat it) and the MAF stage on every cluster, on the same trees and with the same rooting as the solver. The values equal the corresponding entries of optimize(...).meta. Use it to check how hard an instance is before solving.

Parameters:

maf (MAFComputer | None, optional) – Agreement-forest stage. None (default) uses HeuristicMAF.

Returns:

Cluster count and sizes, the MAF size per cluster, and the signature count.

Return type:

MaxMPDStructure

Raises:

  • ForestPDInputError – If the MAF stage does not return an agreement forest (a partition of each cluster’s leaves), or maf is not a MAF computer, or a forest has more than 63 components.

  • ForestPDValueError – From a MAF stage that cannot handle a cluster (e.g. maffe and a polytomy).

Examples

>>> import forestpd as fpd
>>> s = fpd.MaxMPDInstance(["(A:1,B:1,(C:1,D:1):2);", "(A:1,C:1,(B:1,D:1):2);"]).structure()
>>> s.n_clusters, s.maf_sizes, s.max_maf_size
(1, (2,), 2)
evaluate_pd(taxa)[source]#

Objective value of protecting a taxon set.

The expected sum_j w_j PD_{T_j} of the surviving taxa: the given taxa survive surely, every other taxon x independently with probability a(x). Without survival probabilities this is sum_j w_j PD_{T_j}(taxa).

The first call groups the edges of all trees by split once (compiled, linear in the size of the trees); every call is then linear in the size of that split table, independent of the number of trees.

Parameters:

taxa (Iterable[str]) – Taxon labels.

Returns:

The (expected) weighted sum of phylogenetic diversities over all trees.

Return type:

float

Raises:

ForestPDInputError – If a label is unknown, or taxa is a single string rather than a collection.

class forestpd.MaxMPDResult(value, taxa=None, meta=<factory>)[source]#

Result of a MaxMPD optimization.

Variables:

  • value (float) – The optimal objective value: the expected sum_j w_j PD_{T_j} of the surviving taxa when S is protected (sum_j w_j PD_{T_j}(S) without survival probabilities).

  • taxa (set[str] | None) – An optimal taxon-label set, or None when return_taxa was not requested.

  • meta (dict) –

    Diagnostics:

    • "frontier": the optimal value PD(beta) for the budgets beta = 0..B, as float64[min(B, C)+1] with C the total cost of all taxa (a larger budget buys nothing more, so its value is the last entry). The budgets the computed tables do not cover are NaN: see "tables";

    • "tables": which of the manuscript’s tables were computed. "DP" (budget allocated, up to the largest budget) covers every budget up to it; "DP'" (budget saved, up to C minus the smallest budget) covers every budget from the smallest one on; "both" covers the requested budgets from one or the other. The cheapest is chosen: DP’ when the budgets are more than half the total cost;

    • "budget", "n_trees", "n_taxa", "maf" (the MAF stage used);

    • "n_clusters", "cluster_sizes", "maf_sizes", "max_maf_size", "n_signatures": as in MaxMPDStructure (the sizes as lists);

    • timings in seconds: "parse_seconds" (reading and validating the trees, done once when the instance was built, so the same in every result of that instance), "reduction_seconds" (about zero once the instance has cached its reduction), "maf_seconds" and "dp_seconds", whose sum is "solve_seconds".

class forestpd.MaxMPDStructure(n_clusters, cluster_sizes, maf_sizes, max_maf_size, n_signatures, reduction_seconds, maf_seconds)[source]#

The structure that determines the running time, from MaxMPDInstance.structure().

The per-cluster tuples follow the solver’s processing order (children before parents), the order of meta["maf_sizes"] in an optimization result.

Variables:

  • n_clusters (int) – Number of clusters of the common-cluster reduction.

  • cluster_sizes (tuple[int, ...]) – Per cluster, the leaves of its local instance I_C: its own taxa and one pseudo-leaf per child cluster.

  • maf_sizes (tuple[int, ...]) – Per cluster, the size k of its agreement forest.

  • max_maf_size (int) – The largest k: the exponent of the running time O(2^k * n * B^2).

  • n_signatures (int) – sum_C 2^{k_C}, counting the PD_C^+ tables of non-root clusters as well: the number of signatures the DP visits at most (a small budget skips the larger ones).

  • reduction_seconds (float) – Time of the cluster reduction (about zero when the instance had it cached).

  • maf_seconds (float) – Time of the MAF stage.