JaCK & Jill

The content on this page was written by AI under human supervision.

JaCK & Jill is the Python package phyloexact (import phyloexact as px): JaCK for the exact Jukes–Cantor and Kimura computations, Jill for the interval and sampling routines. Given a small alignment, a tree topology and a substitution model, it computes the tree's Bayesian evidence and says which kind of answer it could give: an exact fraction, a proven interval, an estimate with a measured error, or a proof that two trees tie. Every answer carries a machine-readable certificate that a separate routine in the package rechecks against the data.

What it does

Bayesian phylogenetics weighs one tree or one substitution model against another by its evidencemarginal likelihood: the probability of the alignment under a tree and a model, integrated over the branch lengths with a prior; a ratio of two evidences is a Bayes factor. Programs in everyday use estimate the evidence by sampling, without a proven error bound. For four taxa under the Jukes–Cantor model (JC69, all substitutions equally likely) with the package's prior it is a rational number, returned as a Python Fraction:

$$Z(T)=\int_{(0,1)^5} P(D\mid T,x)\,dx_1\cdots dx_5,\qquad x_e=e^{-4t_e/3}.$$

The prior is uniform in each $x_e$, which is an exponential prior of rate $4/3$ on each branch length $t_e$. Under JC69 the likelihood is a polynomial in the $x_e$ with rational coefficients, so $Z$ is rational and a Bayes factor is a ratio of integers. Kimura's two-parameter model (K2P) is exact too on shorter alignments, and protein sequences enter through a four-letter recoding ('JC-SR4').

The package calls the kind of answer its register. R0 is a tie proof: if a permutation of the taxa leaves the site-pattern counts unchanged and maps one topology onto the other, the two evidences are equal under the usual models (JC69, K2P, GTR, mixtures) with any independent branch prior. Nothing is integrated, and every comparison checks for this first. R1 is the exact fraction, with a check that its denominator divides a bound known in advance ($256^N\,\mathrm{lcm}(1,\ldots,N+1)^5$ for $N$ columns under JC69). R2 is an interval with exact rational endpoints proven to contain $Z$, by quadrature in ball arithmetic (needs python-flint); where only the upper end is proven the label is R2-one-sided. R3 is an estimate by importance sampling with scrambled Sobol points, available at any alignment length for JC69, HKY85, GTR with gamma rates and your own log-likelihood (px.PlugIn). Its error figure comes from a calibration record of measured deviations from exact values on small alignments, identified by checksum in the certificate: an indicative accuracy, not a bound. The bundled record was measured on JC69 quartets; HKY85 and GTR estimates reuse it and say so in the certificate, and a model without a usable record, such as general Markov ('GM'), is refused with a message instead of being estimated without one.

register="auto" returns the strongest answer available, based on an import-time probe of the installed libraries and GPU that is recorded in every certificate. Forcing a register that cannot be met raises px.RegisterUnavailable, whose reason states the limit and whose unlock names the smallest change that would satisfy the request; nothing weaker is substituted silently. On a CPU, JC69 is exact for quartets up to 150 usable columns and K2P up to about 35; past 150 columns exact JC69 needs a CUDA GPU with torch, and only for alignments inside a measured size range that the refusal message states. px.verify_certificate rechecks a certificate against the data, recomputing the value where that is cheap and saying what it skipped. It catches mistakes, stale installs and corruption, and is not designed to stop someone who can rewrite both the certificate and the package.

Further calls compare and grade models and check their adequacy (see Routines). px.hill works at genome length. It runs window scans and computes an exact upper bound on the likelihood that any independent-sites model can reach on an alignment of any size. It also checks published likelihood tables against that bound and writes weighted-quartet input files for the species-tree programs wASTRAL and wQFM. px.quintet ranks the 15 rooted five-taxon trees by maximum likelihood with proven bounds, abstaining rather than guessing.

Limits. The package scores a given topology and does not search trees. Exact values and intervals are four-taxon (the quintet engine aside); beyond that it offers tie proofs for up to 32 taxa and the JC69 estimate. HKY85 and GTR have irrational transition probabilities, so no exact value exists and a forced R1 says so. Columns with gaps or ambiguity codes are dropped. The weight files certify inputs to a species-tree program, never the tree.

Examples

Run the self-test suite. We want to know a fresh copy works on this machine before trusting its numbers.

python3 -m phyloexact.validate

The suite prints the capability probe, one line per test with PASS or FAIL and its time, and an overall verdict, and writes VALIDATION_REPORT.json. The tests recompute known exact values on the bundled fixtures, make deliberate edits to certificates and reference files that the verifiers must reject, and confirm that infeasible requests are refused with the right message; a missing optional library is reported as a skip. No GPU or external data is needed, and any failure makes the exit code nonzero.

The smallest exact case, from the suite's JC69 test: one column in which all four taxa carry the same base, then twenty such columns with the exact register forced.

from fractions import Fraction as F
from phyloexact import Dataset, evidence

rep = evidence(Dataset.from_counts({"AAAA": 1}), "12|34", "JC69")
rep.register == "R1" and rep.value == F(103, 4096)        # True
rep = evidence(Dataset.from_counts({"AAAA": 20}), "12|34", "JC69",
               register="R1")
c = rep.certificate
c["denominator_bound"]["divides"]                          # True
len(c["engine"]["sha256"]) == 64                           # True

A constant column has evidence exactly $103/4096$ on 12|34; the same test checks AGCT ($7/4096$) and AACG ($15/4096$). rep.logZ is the floating-point logarithm. In the certificate, denominator_bound records that the denominator divides the bound above, and engine holds the SHA-256 of the exact engine's source, which the verifier compares with the file on disk.

A quartet from a FASTA file, using lines from the package guide's quickstart: one evidence, the tie-proof certificate, a topology comparison, a between-model Bayes factor and a certificate check.

import phyloexact as px
ds  = px.Dataset.from_fasta("quartet.fasta", taxa=["a", "b", "c", "d"])
rep = px.evidence(ds, topology="12|34", model="JC69", register="auto")
rep.value                     # exact Fraction (R1) or float logZ view (R2/R3)
rep.register                  # register actually delivered, e.g. 'R1'
rep.certificate               # machine-checkable dict
px.decidability(ds)           # R0 certificate directly (always available)
px.compare(ds, "13|24", "14|23", model="GTR+G")   # forced tie => R0, no arithmetic
px.compare_models(ds, "12|34", models=("K2P", "JC69"))  # exact between-model BF
px.verify_certificate(rep.certificate, dataset=ds)  # machine-check any cert

taxa= fixes the order, so "12|34" pairs a with b; "13|24" and "14|23" are the other two quartet trees. For a transfer-RNA-length alignment rep.register is 'R1' and rep.value a Fraction with hundreds of digits; on a long one auto gives 'R3' and a floating-point log evidence. px.compare returns a dictionary whose verdict is 'forced-tie' (with log_bf exactly 0 and the symmetry that proves it) or 'computed' (with log_bf, and bf_exact when both sides are exact). px.compare_models returns a ModelComparisonReport whose bf is an exact ratio when both models ran exactly, and px.verify_certificate returns ok and a list of checks marked pass, fail or skipped.

Routines

Data

Evidence and model checking

Certificates and environment

Alignment scale (px.hill)

Five taxa (px.quintet)

Used on this site

Requirements and source

Python 3.10 or later with numpy, scipy, sympy and mpmath; optional python-flint (intervals), numba (compiled sampling kernels), torch with CUDA (large exact JC69). Install with pip install -e . from the root of a jackandjill-dev checkout or put that directory on PYTHONPATH; the package reads checksummed reference files beside it, so it is not built as a wheel, and files under pins/ and phyloexact/quintet/_vendor/ must not be edited. Self-tests: python3 -m phyloexact.validate. Code: BootLoops' jackandjill-dev repository (GitHub organization BootLoops-ai, beside the bootloops-dev toolkit repository rather than inside it), released under the MIT license. The supplementary material on the phylogenetics page is the full manual: every call with its options, returned objects and refusals, the models available for each kind of answer, and cost planning. The package README marks a few interfaces as experimental (the pail_* calls, non-default priors and plug-in models, the genome-scan command-line tools, the quintet tuning options, the estimate route of the adequacy check): they are tested, but their signatures may change in a later release. Authors: Matthew D. Schwartz, Scott V. Edwards, Paul O. Lewis and Claude (Anthropic).

← back to the tools index