#!/usr/bin/env python3 """Phylogenetic quartet (Jukes-Cantor) -- compute the exact Bayesian evidence at runtime from the kinematic inputs (site-pattern counts u, prior shape alpha, topology) and gate it against independent oracles. Self-contained (mpmath + standard library only). NOTHING about the headline results is stored: the exact rational Z(T;u,alpha), the Bayes factor, the GTR probe value, and every agreement digit count are all computed live below. What runs at runtime: Point 1 Z(12|34), Z(13|24), Z(14|23) for the twenty-site synthetic calibration alignment: (i) the 15 JC site-pattern polynomials are built by Felsenstein pruning (exact, fractions); (ii) the product over the alignment is expanded exactly (Kronecker big-integer encoding of the 5-variable coefficient tensor); (iii) each edge is integrated against the Exp-type prior exactly (monomial x^k -> alpha/(k+alpha) in Q), giving Z in Q. Live oracle: independent 5-fold Gauss-Legendre quadrature (exact for the polynomial integrand up to rounding), never touching the expansion route. The exact phylogenetic Bayes factor log[Z(12|34)/Z(13|24)] is formed from the two fresh fractions. Point 2 GTR 3-taxon probe R2: exact rational Z via the tensor-resolvent identity int_0^inf P(t)^{(x)N} lam e^{-lam t} dt = lam (lam I - Q^{(+)N})^{-1} (64x64 exact-rational solve), gated against a live mpmath spectral-eigendecomposition oracle (120 dps). Point 3 +Gamma(k=1) rate heterogeneity: Z = (10-6G)/64 with G = e*E_1(1) the Euler-Gompertz constant; G recomputed live via mp.expint AND via an E_1-free quadrature int_0^inf e^{-t}/(1+t) dt. Retained literals (all held-out oracles, none used in the computation): REF_Z_50D, REF_LOGZ_40D, REF_LOGZ_1324 -- independent packaged decimals from the original run (2026); REF_GTR_SPECTRAL -- spectral oracle value (the packaged spectral reference set; 80 significant digits of the 120-dps run ship with it); REF_G_200D, REF_ZH_200D -- Euler-Gompertz constant and Z_H values (the packaged high-precision reference set; ~101 significant digits of the 200-dps run ship with it, so comparisons saturate near 100 d). Archived full-precision statistics quoted in comments (e.g. the 61.92 d exact-vs-quadrature gate at dps 60+, the 225.4 d GTR two-precision gate) were measured offline in the campaign; the digit counts PRINTED by this script are recomputed here at the working precisions chosen for a short runtime. Interface (evaluate at a DIFFERENT kinematic point / precision, no editing): default run gate demo above, live-gated, plus a dps-doubling check (each live oracle rerun at 2*dps: agreement digits must grow; where a comparison is against a stored oracle string its digit cap is printed). Measured wall time is printed; ~2.5 min at the default --dps 70 on the reference box. --dps D oracle working precision (the doubling check runs at 2*D; the exact rationals are dps-independent). --point 'xxxx:10,xxyy:4,...' [--alpha A] [--topology 12|34] [--dps D] quartet evidence Z(T; u, alpha) at ANY kinematic point: u = counts on the 15 canonical JC pattern labels, alpha = rational prior shape > 0 (finite exact computation, no series radius), T in {12|34, 13|24, 14|23}. Z is EXACT in Q, so any dps is available for free once computed. Cost: the exact tensor has (N+1)^5 slots, N = sum(u) -- measured ~30 s at N=20 (the gate point); keep N <= ~30. The live GL oracle runs when alpha = 1/s (s a positive integer) and costs ((s*N)//2+1)^5 nodes (~12 s at N=20, s=1). --gtr 'pi=2/10,3/10,1/10,4/10;r=1,2,3,1,2,1;lam=1;u=xxx:1,xxy:1,xyz:1' 3-taxon GTR probe at any reversible rate point (pi: 4 positive rationals summing to 1; r: 6 positive rationals; lam > 0; u on labels xxx,xxy,xyx,yxx,xyz with sum(u) <= 3): exact resolvent Z gated against the live spectral oracle. As a library, evaluate(u, alpha, topology, dps) is the callable f(point, dps). """ import argparse import itertools import time from fractions import Fraction from math import gcd import mpmath as mp mp.mp.dps = 220 # default working precision; main() may raise it def digits(a, b): """-log10 |a-b|/|b| : measured agreement in decimal digits.""" a, b = mp.mpf(a), mp.mpf(b) if a == b: return float('inf') return float(-mp.log10(abs((a - b) / b))) def frac_mp(fr): return mp.mpf(fr.numerator) / mp.mpf(fr.denominator) # =========================================================================== # Point 1 machinery: exact Z(T; u, alpha) for the 4-taxon JC quartet # =========================================================================== STATES = (0, 1, 2, 3) def _canonical_label(s): """First-appearance lettering of a state tuple, e.g. (0,0,1,1)->'xxyy'.""" seen, out = {}, [] for si in s: if si not in seen: seen[si] = 'xyzw'[len(seen)] out.append(seen[si]) return ''.join(out) def build_quartet_kernels(): """15 site-pattern polynomials for topology 12|34, by Felsenstein pruning. JC transition along an edge with x = e^{-b} (b = rescaled branch length): P_same = (1+3x)/4, P_diff = (1-x)/4. Returns {label: {(k1..k5) in {0,1}^5: int}} with the dict encoding 256*p_label as an integer-coefficient multilinear polynomial in (x1..x4 pendant, x5 internal), plus the orbit multiplicities. """ def P(i, j): # (constant term, x-coefficient) return (Fraction(1, 4), Fraction(3, 4)) if i == j else \ (Fraction(1, 4), Fraction(-1, 4)) kernels, mults = {}, {} for s in itertools.product(STATES, repeat=4): lab = _canonical_label(s) if lab in kernels: mults[lab] += 1 continue mults[lab] = 1 poly = {} for a in STATES: # internal node adjacent to taxa 1,2 for b in STATES: # internal node adjacent to taxa 3,4 fac = [P(a, s[0]), P(a, s[1]), P(b, s[2]), P(b, s[3]), P(a, b)] for bits in itertools.product((0, 1), repeat=5): c = Fraction(1, 4) # uniform root distribution for e in range(5): c *= fac[e][bits[e]] poly[bits] = poly.get(bits, Fraction(0)) + c k = {} for bits, c in poly.items(): v = c * 256 assert v.denominator == 1 if v: k[bits] = int(v) kernels[lab] = k assert len(kernels) == 15 # partition-of-unity check: sum_k mult_k * 256 p_k == 256 identically tot = {} for lab, k in kernels.items(): for bits, c in k.items(): tot[bits] = tot.get(bits, 0) + mults[lab] * c assert {b: c for b, c in tot.items() if c} == {(0, 0, 0, 0, 0): 256} return kernels, mults KERNELS, PATTERN_MULT = build_quartet_kernels() def permute_dataset(u, tau): """Pattern counts for topology tau*(12|34): read taxon tau[i] at leaf i.""" out = {} # reconstruct a representative state tuple from each label for lab, cnt in u.items(): rep = tuple('xyzw'.index(ch) for ch in lab) s2 = tuple(rep[tau[i] - 1] for i in range(4)) lab2 = _canonical_label(s2) out[lab2] = out.get(lab2, 0) + cnt return out def Z_quartet_exact(u, alpha=Fraction(1)): """Exact Z(12|34; u, alpha) in Q, for any pattern counts u and any rational prior shape alpha (edge prior on x: alpha x^{alpha-1} dx on [0,1], i.e. Exp(alpha) on the rescaled branch length b = -log x; monomial integral: int_0^1 x^k alpha x^{alpha-1} dx = alpha/(k+alpha)). The product of the N site polynomials is expanded EXACTLY via a Kronecker big-integer encoding: the 5-variable coefficient tensor (degree <= N per axis) is packed into one Python integer, one slot of `bslot` bits per coefficient, so polynomial multiplication becomes shifted big-integer addition. Slot width is set from the product of the kernels' L1 norms, which bounds every intermediate coefficient. """ N = sum(u.values()) nax = N + 1 # slots per axis strides = [nax ** i for i in range(5)] bound = 1 for lab, e in u.items(): bound *= sum(abs(c) for c in KERNELS[lab].values()) ** e bslot = ((bound.bit_length() + 2 + 7) // 8) * 8 # byte-aligned slot # ---- exact product over the alignment (Kronecker encoding) ---- V = 1 for lab, e in u.items(): terms = [(sum(bits[i] * strides[i] for i in range(5)) * bslot, c) for bits, c in KERNELS[lab].items()] for _ in range(e): V = sum((V << sh) * c for sh, c in terms) # ---- decode the coefficient tensor (signed slots via offset trick) ---- nslots = nax ** 5 half = 1 << (bslot - 1) offset = half * ((1 << (bslot * nslots)) - 1) // ((1 << bslot) - 1) raw = (V + offset).to_bytes(nslots * bslot // 8, 'little') bb = bslot // 8 coeffs = [int.from_bytes(raw[i * bb:(i + 1) * bb], 'little') - half for i in range(nslots)] # ---- exact edge integrals: x^k -> alpha/(k+alpha) ---- w = [Fraction(alpha) / (k + Fraction(alpha)) for k in range(nax)] L = 1 for wk in w: L = L * wk.denominator // gcd(L, wk.denominator) wi = [int(wk * L) for wk in w] # integer-scaled weights cur = coeffs for _axis in range(5): # contract one axis at a time cur = [sum(cur[j + k] * wi[k] for k in range(nax) if cur[j + k]) for j in range(0, len(cur), nax)] return Fraction(cur[0], 256 ** N * L ** 5) # ---- live oracle: independent 5-fold Gauss-Legendre quadrature ----------- def _gl_nodes(m, dps): """Gauss-Legendre nodes/weights on [0,1] (Newton on Legendre P_m).""" mp.mp.dps = dps xs, ws = [], [] for k in range(1, m + 1): x = mp.cos(mp.pi * (k - mp.mpf(1) / 4) / (m + mp.mpf(1) / 2)) for _ in range(100): p0, p1 = mp.mpf(1), x for j in range(2, m + 1): p0, p1 = p1, ((2 * j - 1) * x * p1 - (j - 1) * p0) / j dp = m * (x * p1 - p0) / (x * x - 1) dx = -p1 / dp x += dx if abs(dx) < mp.mpf(10) ** (-dps + 2): break xs.append((x + 1) / 2) ws.append(1 / ((1 - x * x) * dp * dp)) return xs, ws def Z_quartet_gl(u, dps=70, s=1): """Live quadrature oracle for alpha = 1/s, s integer >= 1 (s=1: uniform in x): after the substitution x = y^s, int_0^1 f(x) alpha x^{alpha-1} dx = int_0^1 f(y^s) dy (alpha = 1/s), the integrand is a polynomial of per-axis degree s*N, so an m = (s*N)//2+1 point Gauss-Legendre rule is EXACT up to rounding. Independent of the expansion route: never forms the coefficient tensor, just evaluates the pruning likelihood at m^5 nodes. Uses the pruning factorization p_s = (1/4) [ (1-x5)/4 * (sum_a A_a)(sum_b B_b) + x5 * sum_a A_a B_a ], A_a = P(a,s1;x1) P(a,s2;x2), B_b = P(b,s3;x3) P(b,s4;x4). """ dps_save = mp.mp.dps N = sum(u.values()) m = (s * N) // 2 + 1 xs, ws = _gl_nodes(m, dps) mp.mp.dps = dps if s != 1: xs = [y ** s for y in xs] # x = y^s; weights unchanged (dy measure) pats = [(tuple('xyzw'.index(ch) for ch in lab), e) for lab, e in u.items()] def P(i, j, x): return (1 + 3 * x) / 4 if i == j else (1 - x) / 4 quarter = mp.mpf(1) / 4 S = mp.mpf(0) for i1 in range(m): x1 = xs[i1] for i2 in range(m): x2 = xs[i2] w12 = ws[i1] * ws[i2] Avec = {} for (s, e) in pats: key = (s[0], s[1]) if key not in Avec: Avec[key] = [P(a, key[0], x1) * P(a, key[1], x2) for a in range(4)] for i3 in range(m): x3 = xs[i3] for i4 in range(m): x4 = xs[i4] w1234 = w12 * ws[i3] * ws[i4] Bvec = {} for (s, e) in pats: key = (s[2], s[3]) if key not in Bvec: Bvec[key] = [P(b, key[0], x3) * P(b, key[1], x4) for b in range(4)] pre = [] for (s, e) in pats: A = Avec[(s[0], s[1])] B = Bvec[(s[2], s[3])] pre.append((sum(A) * sum(B), sum(A[i] * B[i] for i in range(4)), e)) for i5 in range(m): x5 = xs[i5] q = (1 - x5) / 4 val = w1234 * ws[i5] for (SASB, dAB, e) in pre: val *= (quarter * (q * SASB + x5 * dAB)) ** e S += val mp.mp.dps = dps_save return S # =========================================================================== # Point 1: quartet Z, twenty-site synthetic calibration alignment, alpha=1 # =========================================================================== U_REF = {'xxxx': 10, 'xxyy': 4, 'xyxy': 2, 'xyyx': 2, 'xxxy': 1, 'xyzw': 1} # held-out packaged oracles (independent original run, 2026). NOT used above. REF_Z_50D = "8.5983925891730766062850061830326329145772650106521e-39" REF_LOGZ_40D = "-87.64924334921352951083970264655772400877" REF_LOGZ_1324 = "-90.83748304235128116508526548224139936558" # =========================================================================== # Point 2: GTR probe R2 (3-taxon, generic rates) -- tensor-resolvent, exact # =========================================================================== IDX = {(0, 1): 0, (0, 2): 1, (0, 3): 2, (1, 2): 3, (1, 3): 4, (2, 3): 5} REPS3 = {'xxx': (0, 0, 0), 'xxy': (0, 0, 1), 'xyx': (0, 1, 0), 'yxx': (1, 0, 0), 'xyz': (0, 1, 2)} def Q_gtr(pi, r): """GTR generator: Q_ij = pi_j r_{ij} (i != j), rows sum to zero.""" Q = [[Fraction(0)] * 4 for _ in range(4)] for i in range(4): for j in range(4): if i != j: Q[i][j] = pi[j] * r[IDX[(min(i, j), max(i, j))]] for i in range(4): Q[i][i] = -sum(Q[i][j] for j in range(4) if j != i) return Q def _kron_sum(Q, N): """Q^{(+)N} = sum_k I x..x Q x..x I as a 4^N x 4^N Fraction matrix.""" dim = 4 ** N M = [[Fraction(0)] * dim for _ in range(dim)] for I in range(dim): iv = [(I // 4 ** k) % 4 for k in range(N)] for k in range(N): for jp in range(4): if jp == iv[k]: M[I][I] += Q[iv[k]][iv[k]] else: M[I][I + (jp - iv[k]) * 4 ** k] += Q[iv[k]][jp] return M def _solve_exact(A, b): """Gaussian elimination over Fraction.""" n = len(A) M = [row[:] + [b[i]] for i, row in enumerate(A)] for col in range(n): piv = next(r for r in range(col, n) if M[r][col] != 0) M[col], M[piv] = M[piv], M[col] inv = Fraction(1) / M[col][col] M[col] = [x * inv for x in M[col]] for r in range(n): if r != col and M[r][col] != 0: f = M[r][col] M[r] = [M[r][j] - f * M[col][j] for j in range(n + 1)] return [M[i][n] for i in range(n)] def Z_gtr_resolvent(pi, r, u, lam): """EXACT rational Z for the 3-taxon star tree under GTR + Exp(lam) prior, via int_0^inf P(t)^{(x)N} lam e^{-lam t} dt = lam (lam I - Q^{(+)N})^{-1} -- a 4^N x 4^N exact-rational linear solve per distinct edge column.""" pi = [Fraction(p) for p in pi] r = [Fraction(x) for x in r] lam = Fraction(lam) Qf = Q_gtr(pi, r) sites = [] for l, n in u.items(): sites.extend([REPS3[l]] * n) Ns = len(sites) dim = 4 ** Ns QN = _kron_sum(Qf, Ns) cols, Re = {}, [] for e in range(3): tgt = sum(sites[i][e] * 4 ** i for i in range(Ns)) if tgt not in cols: A = [[(-QN[i][j] if i != j else lam - QN[i][i]) for j in range(dim)] for i in range(dim)] b = [Fraction(0)] * dim b[tgt] = lam cols[tgt] = _solve_exact(A, b) Re.append(cols[tgt]) Z = Fraction(0) for a in range(dim): av = [(a // 4 ** k) % 4 for k in range(Ns)] w = Fraction(1) for i in range(Ns): w *= pi[av[i]] for e in range(3): w *= Re[e][a] Z += w return Z def Z_gtr_spectral(pi, r, u, lam, dps): """Live independent oracle: mpmath symmetric eigendecomposition of D Q D^{-1} and spectral sum over internal/leaf state assignments.""" dps_save = mp.mp.dps mp.mp.dps = dps def fm(x): f = Fraction(x) return mp.mpf(f.numerator) / mp.mpf(f.denominator) pim = [fm(p) for p in pi] rm = [fm(x) for x in r] lam = fm(lam) Q = mp.zeros(4, 4) for i in range(4): for j in range(4): if i != j: Q[i, j] = pim[j] * rm[IDX[(min(i, j), max(i, j))]] for i in range(4): Q[i, i] = -sum(Q[i, j] for j in range(4) if j != i) D = mp.diag([mp.sqrt(p) for p in pim]) Di = mp.diag([1 / mp.sqrt(p) for p in pim]) ev, U = mp.eigsy(D * Q * Di) A = Di * U B = U.T * D sites = [] for l, n in u.items(): sites.extend([REPS3[l]] * n) Ns = len(sites) Z = mp.mpf(0) for av in itertools.product(range(4), repeat=Ns): w = mp.mpf(1) for i in range(Ns): w *= pim[av[i]] for e in range(3): Se = mp.mpf(0) for lv in itertools.product(range(4), repeat=Ns): c = mp.mpf(1) es = mp.mpf(0) for i in range(Ns): c *= A[av[i], lv[i]] * B[lv[i], sites[i][e]] es += ev[lv[i]] Se += c * lam / (lam - es) w *= Se Z += w mp.mp.dps = dps_save return Z # R2 probe kinematics (generic rate point, eigenvalue cubic irreducible/Q): PI_R2 = ['2/10', '3/10', '1/10', '4/10'] R_R2 = ['1', '2', '3', '1', '2', '1'] U_R2 = {'xxx': 1, 'xxy': 1, 'xyz': 1} # held-out packaged spectral value (spectral reference set, 120-dps original # run; 80 significant digits ship here) REF_GTR_SPECTRAL = ("0.000000588272335916668608716289785653467205548337008110" "41993035628155738563417962236587") # =========================================================================== # Point 3: +Gamma(k=1) rate heterogeneity, Z = (10 - 6G)/64 # =========================================================================== # held-out packaged values (high-precision reference set, 200-dps original # run; ~101 significant digits ship here, so agreements saturate near 100 d) REF_G_200D = ("0.59634736232319407434107849936927937607417786015254878157" "348491048232721911487441747043049709361276034") REF_ZH_200D = ("0.10034243478220055553052389068413005849304582561069855172" "748578964228182320798052336214714089747380371") # =========================================================================== # Evaluation interface: f(kinematic point, dps) + gate-demo driver # =========================================================================== TOPOLOGIES = {'12|34': (1, 2, 3, 4), '13|24': (1, 3, 2, 4), '14|23': (1, 4, 2, 3)} DEMO_DPS = 70 # default oracle precision for the gate demo def _cap(sref): """Significant digits carried by a stored oracle string: agreement against it saturates here; the live-vs-live gates have no such cap.""" return len(sref.split('e')[0].replace('-', '').replace('.', '').lstrip('0')) def evaluate(u, alpha=Fraction(1), topology='12|34', dps=50, oracle=True): """f(kinematic point, dps) for the quartet evidence Z(topology; u, alpha). Domain: u = site-pattern counts (nonnegative ints on the 15 canonical JC labels = KERNELS keys), N = sum(u) >= 1; alpha = rational prior shape > 0 (the computation is a finite exact sum -- no series radius); topology in TOPOLOGIES. Z is EXACT in Q; dps sets only the returned mpf precision and the live-oracle precision. Cost: (N+1)^5 tensor slots, measured ~30 s at N=20; live GL oracle needs alpha = 1/s, s a positive integer. """ if topology not in TOPOLOGIES: raise ValueError("topology must be one of %s" % sorted(TOPOLOGIES)) bad = [l for l in u if l not in KERNELS] if bad or not u or any(int(c) != c or c < 0 for c in u.values()): raise ValueError("u must map canonical labels %s to nonnegative ints" " (bad labels: %s)" % (sorted(KERNELS), bad)) alpha = Fraction(alpha) if alpha <= 0: raise ValueError("alpha must be a positive rational") mp.mp.dps = max(mp.mp.dps, dps + 20) uu = permute_dataset(u, TOPOLOGIES[topology]) t0 = time.time() Zf = Z_quartet_exact(uu, alpha) res = {'Z_frac': Zf, 'Z': frac_mp(Zf), 'logZ': mp.log(frac_mp(Zf)), 't_exact_s': time.time() - t0, 'oracle': None} if oracle and (1 / alpha).denominator == 1: s = int(1 / alpha) t0 = time.time() Zgl = Z_quartet_gl(uu, dps=dps, s=s) res['oracle'] = {'Z_gl': Zgl, 'digits': digits(Zgl, res['Z']), 'dps': dps, 't_gl_s': time.time() - t0} return res def demo(dps=DEMO_DPS): """Default run: the three gate points computed live, gated against the live oracles and the held-out packaged strings, plus the dps-doubling check (each live oracle rerun at 2*dps; agreement digits must grow).""" dps2 = 2 * dps mp.mp.dps = max(220, dps2 + 40) T0 = time.time() print("Gate demo: live-oracle dps=%d, dps-doubling check at dps=%d" % (dps, dps2)) print("\n== Point 1: quartet Z, twenty-site synthetic calibration" " alignment, alpha=1 ==") Zs = {} for topo in ('12|34', '13|24', '14|23'): t0 = time.time() Zs[topo] = Z_quartet_exact(permute_dataset(U_REF, TOPOLOGIES[topo])) print(" Z(%s) exact rational computed live: %.1f s" % (topo, time.time() - t0)) z = frac_mp(Zs['12|34']) print(" Z(12|34) =", mp.nstr(z, 50)) print(" Z (packaged reference) =", REF_Z_50D) print(" agreement : %.1f d (stored string caps at %d d)" % (digits(z, mp.mpf(REF_Z_50D)), _cap(REF_Z_50D))) lz, lz2 = mp.log(z), mp.log(frac_mp(Zs['13|24'])) print(" log Z(12|34) =", mp.nstr(lz, 40)) print(" log Z (packaged, 40 d) =", REF_LOGZ_40D) print(" agreement : %.1f d (stored string caps at %d d)" % (digits(lz, mp.mpf(REF_LOGZ_40D)), _cap(REF_LOGZ_40D))) print(" log Z(13|24) =", mp.nstr(lz2, 40)) print(" log Z(13|24) (packaged, 40 d) =", REF_LOGZ_1324) print(" agreement : %.1f d" % digits(lz2, mp.mpf(REF_LOGZ_1324))) print(" Z(14|23) == Z(13|24) exactly :", Zs['14|23'] == Zs['13|24'], "(dataset symmetry; both fractions computed independently)") print(" Bayes factor log[Z(12|34)/Z(13|24)] =", mp.nstr(lz - lz2, 40)) print(" live GL quadrature oracle (independent route), dps-doubling:") for dd in (dps, dps2): t0 = time.time() S = Z_quartet_gl(U_REF, dps=dd) print(" dps=%4d: agreement vs exact = %7.2f d (%.1f s)" % (dd, digits(S, z), time.time() - t0)) print(" (live-vs-live gate: no stored-float cap; digits track dps)") print("\n== Point 2: GTR 3-taxon probe R2 (tensor-resolvent, exact) ==") t0 = time.time() Zg = Z_gtr_resolvent(PI_R2, R_R2, U_R2, 1) zg = frac_mp(Zg) print(" Z exact rational (%.1f s) =" % (time.time() - t0), mp.nstr(zg, 50)) print(" Z (spectral oracle, packaged) =", REF_GTR_SPECTRAL) print(" agreement : %.1f d (stored string caps at %d d)" % (digits(zg, mp.mpf(REF_GTR_SPECTRAL)), _cap(REF_GTR_SPECTRAL))) print(" live spectral oracle, dps-doubling:") for dd in (dps, dps2): t0 = time.time() Zsp = Z_gtr_spectral(PI_R2, R_R2, U_R2, 1, dd) print(" dps=%4d: agreement vs exact = %7.2f d (%.1f s)" % (dd, digits(Zsp, zg), time.time() - t0)) print("\n== Point 3: +Gamma(k=1) 3-taxon, Z = (10 - 6G)/64, G = e*E_1(1) ==") G = mp.e * mp.expint(1, 1) # live, at ambient (>= 2*dps+40) dps Zh = (10 - 6 * G) / 64 print(" G (mpmath e*E_1(1)) =", mp.nstr(G, 50)) print(" G (packaged, 200-dps run) =", REF_G_200D[:52] + "...") print(" agreement on G : %.1f d (stored string caps at %d d)" % (digits(G, mp.mpf(REF_G_200D)), _cap(REF_G_200D))) print(" Z (this run) =", mp.nstr(Zh, 50)) print(" agreement on Z vs packaged : %.1f d (stored string caps at %d d)" % (digits(Zh, mp.mpf(REF_ZH_200D)), _cap(REF_ZH_200D))) print(" live E_1-free quadrature int_0^inf e^-t/(1+t) dt, dps-doubling:") for dd in (dps, dps2): t0 = time.time() with mp.workdps(dd): Gq = mp.quad(lambda t: mp.exp(-t) / (1 + t), [0, mp.inf]) print(" dps=%4d: agreement vs e*E_1(1) = %7.2f d (%.1f s)" % (dd, digits(Gq, G), time.time() - t0)) print("\nTotal wall time: %.1f s" % (time.time() - T0)) def _parse_u(text, allowed): u = {} for item in text.split(','): lab, _, cnt = item.strip().partition(':') if lab not in allowed: raise SystemExit("unknown pattern label %r (allowed: %s)" % (lab, ','.join(sorted(allowed)))) u[lab] = u.get(lab, 0) + int(cnt) return u def run_point(args): u = _parse_u(args.point, KERNELS) if args.point else dict(U_REF) alpha = Fraction(args.alpha or '1') dps = args.dps or 50 T0 = time.time() print("Quartet point: u=%s (N=%d), alpha=%s, topology=%s, dps=%d" % (u, sum(u.values()), alpha, args.topology, dps)) res = evaluate(u, alpha, args.topology, dps=dps, oracle=not args.no_oracle) Zf = res['Z_frac'] if len(str(Zf.numerator)) + len(str(Zf.denominator)) <= 200: print(" Z exact =", Zf) else: print(" Z exact = (%d-digit numerator)/(%d-digit denominator)" % (len(str(Zf.numerator)), len(str(Zf.denominator)))) print(" Z =", mp.nstr(res['Z'], dps)) print(" log Z =", mp.nstr(res['logZ'], dps)) print(" exact route wall time: %.1f s" % res['t_exact_s']) o = res['oracle'] if o: print(" live GL oracle (dps=%d): agreement = %.2f d (%.1f s)" % (o['dps'], o['digits'], o['t_gl_s'])) elif not args.no_oracle: print(" (live GL oracle skipped: it is exact only for alpha = 1/s," " s a positive integer)") print("Total wall time: %.1f s" % (time.time() - T0)) def run_gtr(spec, dps): fields = {} for part in spec.split(';'): k, _, v = part.strip().partition('=') fields[k.strip()] = v.strip() pi = [Fraction(x) for x in fields.get('pi', '1/4,1/4,1/4,1/4').split(',')] r = [Fraction(x) for x in fields.get('r', '1,1,1,1,1,1').split(',')] lam = Fraction(fields.get('lam', '1')) u = _parse_u(fields['u'], REPS3) if 'u' in fields else dict(U_R2) if len(pi) != 4 or sum(pi) != 1 or min(pi) <= 0: raise SystemExit("pi must be 4 positive rationals summing to 1") if len(r) != 6 or min(r) <= 0 or lam <= 0: raise SystemExit("r must be 6 positive rationals and lam > 0") if not 1 <= sum(u.values()) <= 3: raise SystemExit("sum(u) must be 1..3 (the exact 4^N resolvent solve" " is measured at ~2.5 s for N=3 and grows as 4^{3N})") T0 = time.time() print("GTR 3-taxon point: pi=%s r=%s lam=%s u=%s, dps=%d" % ([str(p) for p in pi], [str(x) for x in r], lam, u, dps)) mp.mp.dps = max(mp.mp.dps, dps + 40) t0 = time.time() Zg = Z_gtr_resolvent(pi, r, u, lam) print(" Z exact = %d / %d (%.1f s)" % (Zg.numerator, Zg.denominator, time.time() - t0)) zg = frac_mp(Zg) print(" Z =", mp.nstr(zg, min(dps, 50))) t0 = time.time() Zsp = Z_gtr_spectral(pi, r, u, lam, dps) print(" live spectral oracle (dps=%d): agreement = %.2f d (%.1f s)" % (dps, digits(Zsp, zg), time.time() - t0)) print("Total wall time: %.1f s" % (time.time() - T0)) def main(): ap = argparse.ArgumentParser( description="Exact phylogenetic Bayesian evidence: gate demo (default)" " + arbitrary-point evaluator (see module docstring for" " domain limits and measured costs).") ap.add_argument('--point', metavar='U', help="quartet site-pattern counts, e.g." " 'xxxx:10,xxyy:4,xyxy:2,xyyx:2,xxxy:1,xyzw:1'") ap.add_argument('--alpha', metavar='A', help='rational prior shape > 0 (default 1)') ap.add_argument('--topology', default='12|34', choices=sorted(TOPOLOGIES), help='quartet topology (default 12|34)') ap.add_argument('--dps', type=int, help='oracle/output precision (demo default %d; point' ' mode default 50)' % DEMO_DPS) ap.add_argument('--no-oracle', action='store_true', help='point mode: skip the live GL oracle') ap.add_argument('--gtr', metavar='SPEC', help="GTR probe point, e.g. 'pi=2/10,3/10,1/10,4/10;" "r=1,2,3,1,2,1;lam=1;u=xxx:1,xxy:1,xyz:1'") args = ap.parse_args() if args.gtr: run_gtr(args.gtr, args.dps or 60) elif args.point or args.alpha or args.topology != '12|34' or args.no_oracle: run_point(args) else: demo(dps=args.dps or DEMO_DPS) if __name__ == '__main__': main()