Wayfinder

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

Wayfinder integrates the system of differential equations that a family of Feynman integrals satisfies, carrying high-precision values from a point where they are known to the point where they are wanted, and then extracts their expansion in the regulator $\varepsilon$. A second part, the connection sampler, works the other way: it reconstructs the constant matrices of the differential equation from numerical values of the integrals. Two smaller modules cover neighboring jobs: epslimit takes values of an integral that another program has already computed at several $\varepsilon$ and extracts the $\varepsilon\to0$ expansion from them, and flintexport does exact rational-function algebra in bulk with the FLINT library and writes the exact system files that the loaders read. Input is a saved system (or a small Python object describing one) plus boundary values. Output is the transported vector with a measured truncation error, $\varepsilon$-expansion coefficients with a measured number of reliable digits, and pass/fail comparisons against stored reference values.

What it does

The master integrals of a family obey a linear system $y'(x)=A(x,\varepsilon)\,y(x)$ in a kinematic variable $x$, where $\varepsilon$ is the dimensional regularizationspace-time is continued to d = 4 − 2ε dimensions; results are Laurent series in ε parameter and $A$ is a matrix of rational functions. If the integrals are known at one $x$, for instance a limit where they reduce to vacuum integrals, integrating the system gives them at any other. Wayfinder does this at a fixed numerical $\varepsilon$ with high-order Taylor steps. The step size follows the distance to the nearest singular point of $A$. A step is accepted only when a bound on the neglected Taylor tail is below the requested precision, and the worst such bound over the path is returned as trunc_worst. Complex waypoints take the path around singular points on the real axis; a straight path that runs into one stops with an error asking for a detour. The number of digits is an explicit argument everywhere, and the transport, Frobenius, boundary and $\varepsilon$-expansion routines never change the global mpmath precision (the epslimit module and the connection sampler do set it, as noted below).

A Taylor march cannot end on a regular-singular pointa point where A has a simple pole; solutions there behave like powers (x−x₀)^λ times polynomials in log(x−x₀), so there Wayfinder builds the local Frobenius basis (power-and-logarithm series, including exponents that differ by integers). It then matches the transported vector onto that basis by least squares and returns the residual and a condition-number estimate. A scalar variant handles one higher-order equation given as an exact recurrence, and a two-branch variant handles apparent higher poles. Frobenius matching requires real $\varepsilon$; plain transport accepts complex $\varepsilon$.

The $\varepsilon$-expansion is recovered by repeating the fixed-$\varepsilon$ computation at several $\varepsilon$ and solving for the Laurent coefficients. The nodes sit either on a small circle around $\varepsilon=0$ (default radius $10^{-3}$, 32 nodes) or on a geometric grid of small real values with extra nodes kept back for checking. Every extraction is done twice and the number of digits on which the two agree is returned (min_agree_digits, self_digits); quote that number, not the working precision. Certified variants take the radius from the exact $\varepsilon$-singularities of the system and raise EpsFanCertifyError when the bound cannot be met.

Around this core sit loaders for three on-disk formats, boundary values in AMFlow's normalization, and a two-precision comparison against stored reference digits, in which a value passes only if it matches and does not move when the precision is raised. Smaller modules provide checksum manifests of input files, an adaptive quadrature routine that refines until successive levels agree and otherwise raises, and a certified tail bound for truncated Bessel-moment integrals. A separate module builds the exact connection for one-parameter integrals over a conic $Y^2=q(x,t)$ and refuses higher-degree curves with ConicScopeError.

The connection sampler (sampler.py) covers the case where the equation itself is missing. In a canonical basisa choice of master integrals in which ε factors out of the differential equation and only logarithmic one-forms appear the system reads

$$ dJ \;=\; \varepsilon \Big(\sum_a A_a\, d\log a\Big)\, J, $$

with constant rational matrices $A_a$, one per letter $a$. The sampler evaluates the integrals and their derivatives with AMFlow at a few exact rational kinematic points, by central finite differences or with AMFlow's exact derivative mode. It then solves the linear equations these samples impose on the entries of the $A_a$, rounds each entry to a fraction with PSLQ, and checks the result at points not used in the fit. The output connection.json holds the matrices in sparse form plus diagnostics (points used, condition number, fit residual digits, minimum digits at the check points, entries that failed to rationalize); Ansatzer reads it through --connection.

The epslimit module is for the case where the values already exist. Another program, AMFlow, pySecDec or anything else, has evaluated an integral at a grid of $\varepsilon_i$ of its own choosing, and what is wanted are the coefficients $c_{k,j}$ of $f(\varepsilon)=\sum_{k\ge k_{\min}}\sum_{j\ge0}c_{k,j}\,\varepsilon^k\log^j\varepsilon$, usually the finite part $c_{0,0}$. Samples may be given as numbers, decimal strings or Arb-style [mid +/- rad] strings, and load_amflow_grid collects them from AMFlow result files, skipping any the engine marked as failed. The pole order $k_{\min}$ is read from the slope of $\log|f|$ against $\log|\varepsilon|$. A power of $\log\varepsilon$ is accepted only when refitting with each sample left out in turn improves the fit decisively, which separates a genuine logarithm from round-off in the input, and the working precision is set from the spread and clustering of the nodes. Up to four methods run as the layout of the grid allows: a rescaled Vandermonde solve, Richardson elimination on a geometric ladder of real nodes, rational extrapolation, and a discrete Fourier transform when the nodes lie on a circle around $\varepsilon=0$. The reported achieved_digits is the number of digits on which these methods agree with one another. What can be recovered is fixed by the grid rather than by the algorithm: three samples give six to eight digits of the finite part however they are processed, and thirty digits or more call for six to eight samples on a geometric ladder of ratio 2 to 4, computed at two to three times the target precision. A second routine, richardson_boundary, sums a series at its radius of convergence from a list of partial sums when the powers appearing in the tail are known (half-integers by default, the square-root case). All of this differs from the $\varepsilon$-expansion above, where Wayfinder evaluates the system at nodes it chooses; here the grid is taken as given, and the two share no code. extrapolate sets the global mpmath precision for the duration of the call and restores it before returning; the connection sampler's routines also set it, to the dps they are given, and leave it there.

The flintexport module belongs to the step before transport, when a differential equation or a large set of reduction rules is being assembled symbolically and SymPy's cancel and together become the bottleneck. It walks each SymPy expression once, holding every intermediate value as a pair of exact multivariate polynomials over the rationals in FLINT (through python-flint), removing common factors on addition and cross-canceling on multiplication. Substitutions such as $d\to4-2\varepsilon$ or fixed kinematic values are applied at the leaves. canonical_terms then returns the numerator and denominator as sorted lists of terms with coprime integer coefficients and a positive leading denominator coefficient, a form that can be written to a file or compared byte for byte; this is the format load_monomial_json reads. Any number of variables is allowed. Exponents must be integers, and a fractional or symbolic power raises ValueError instead of being rounded. The speed comes from never leaving FLINT between the leaves and the export, so mixing SymPy arithmetic into the middle of a walk loses it. Counterweight carries the Julia counterpart, MpolyFeed.jl.

Wayfinder does not do integration-by-parts reduction or find a canonical basis; see Kira and Counterweight. Transport works at one $\varepsilon$ at a time by design. The optional python-flint step kernel (backend="acb") is many times faster on systems with hundreds of masters but, like the default, reports a truncation bound rather than interval enclosures. The kira_target.m loader covers a subset of Kira's grammar; unsupported constructs there, and unsupported cases anywhere in the package, raise NotImplementedError naming what is missing instead of returning a guess.

Examples

Run the end-to-end control.

python3 tests/control_hypergeom.py     # end-to-end positive control, ~6 min
python3 tests/control_tadpole.py       # boundary vs direct Gamma, seconds

The first script builds the $2\times2$ system of ${}_2F_1(\varepsilon,-2\varepsilon;1-3\varepsilon;x)$ inline and transports it at $\varepsilon=1/10$ from $x=1/2$ to $x=-1$, around the singular point at $0$ through the waypoint $-0.6i$, at 100 and at 140 digits. The result is compared with 90 stored digits of mpmath.hyp2f1 at $x=-1$, a value the transport never used. It then extracts the $\varepsilon$-coefficients $k=-2\ldots6$ and compares them with an independent Taylor expansion (the second-order one is exactly $\pi^2/6$). Each stage prints its measured digits; exit code 0 means every bar was met. The second script checks tadpole against the Gamma-function closed form at 150 digits.

Transport your own system from Python. Any object with .n, .var, .A(x, eps, dps) and .meta is accepted; this is the control's hypergeometric system with its singular points declared.

from mpmath import mp, mpc
from wayfinder import transport_fixed_eps

class Hyp2F1System:
    n = 2
    var = "x"
    meta = {"source": "inline", "singular_points": [0, 1]}
    def A(self, x, eps, dps):
        with mp.workdps(dps):
            xv, e = mpc(x), mpc(eps)
            den = xv * (1 - xv)
            return [[mpc(0), mpc(1)],
                    [-2 * e * e / den, ((1 - e) * xv - (1 - 3 * e)) / den]]

with mp.workdps(130):                              # boundary vector (w, w') at x = 1/2
    a, b, c = mpc("0.1"), mpc("-0.2"), mpc("0.7")
    y0 = [mp.hyp2f1(a, b, c, 0.5), a * b / c * mp.hyp2f1(a + 1, b + 1, c + 1, 0.5)]

y, diag = transport_fixed_eps(Hyp2F1System(), "0.1", "0.5", "-1", y0, 100,
                              path=[mpc(0, "-0.6")], return_diag=True)

y is $(w,w')$ at $x=-1$ as mpc numbers to the requested 100 digits; diag carries steps, trunc_worst, trunc_worst_log10 and min_hfrac, to be stored with the result. For a system loaded from disk, call desys.enable_fast_path(eps, wp=250) first: it attaches exact Taylor coefficients of $A$ and the singular points, which is much faster on large matrices.

Fit a canonical connection from samples. The one-loop massless box target comes with the package in targets/; this needs an amflow-cpp build (AMFLOW_CLI) and the Gatekeeper and Counterweight packages as siblings under tools/.

python sampler.py targets/box1l_conn.json \
    --n-fit 3 --n-heldout 2 --goal-digits 40 --eps-power 3 \
    --parallel 6 --out out/box1l/connection.json

Three fit points and two verification points cost 25 AMFlow calls, about two minutes at --parallel 6; a cached rerun takes a fraction of a second. In the basis $(\varepsilon(1-2\varepsilon)\,\mathrm{bub}_s,\ \varepsilon(1-2\varepsilon)\,\mathrm{bub}_t,\ \varepsilon^2 st\,\mathrm{box})$ the recovered matrices are

A[s]   = [[-1,0,0],[0,0,0],[2,0,-1]]
A[t]   = [[0,0,0],[0,-1,0],[0,2,-1]]
A[s+t] = [[0,0,0],[0,0,0],[-2,-2,1]]

with a 31-digit fit residual, 32 digits at the verification points, nothing that failed to rationalize, and $\sum_a A_a=-1$ as the scaling check. From Python: from wayfinder.sampler import connection_once; conn = connection_once("targets/box1l_conn.json").

Take the ε → 0 limit of values computed elsewhere. Three small files distributed with the tests hold one sample each of $f(\varepsilon)=7/2-3\varepsilon+11\varepsilon^2$, at $\varepsilon=1/1000$, $1/4000$ and $1/8000$, in the layout of an AMFlow result file.

from wayfinder.epslimit import extrapolate, load_amflow_grid

paths = [f"tests/fixtures_epslimit/series_quad_eps_1_{n}.json" for n in (1000, 4000, 8000)]
samples = load_amflow_grid(paths, part="re")     # [(eps, value), ...], one integral
r = extrapolate(samples, max_log=0)
print(r.finite_part, r.c(1), r.c(2))             # 7/2, -3, 11
print(r.k_min, r.log_order, r.achieved_digits)

r is an ExtrapResult (fields under Routines). Because this $f$ is exactly quadratic, three samples are enough and the test requires $c_0$, $c_1$ and $c_2$ to 40, 35 and 30 digits; for a general function three samples would give six to eight. max_log=0 says no logarithms are expected; the default of 3 lets the routine test for them. The fourth fixture, series_pole_ladder_8pt.json, is an eight-point ladder of $\zeta(3)/\varepsilon+\pi^2+\varepsilon/(1+\varepsilon)$ with one failed sample that the loader must drop; there r.k_min comes back $-1$. From the shell, python3 tools/wayfinder/epslimit.py FILE ... --part re prints the coefficients and ends with a line beginning FINITE PART c_0 =.

Exact rational-function algebra with flintexport. This is the pattern from the module's self-test: SymPy symbols are mapped into the polynomial ring $\mathbb{Q}[\varepsilon,s_2,s_{12}]$, with $d$ sent to $4-2\varepsilon$, and an expression is evaluated there exactly (python-flint required).

import sympy as sp
from wayfinder.flintexport import mk_ctx, gen, poly, make_env, walk, canonical_terms

d, s2, s12 = sp.symbols("d s2 s12")
ctx = mk_ctx(("eps", "s2", "s12"))                      # variables of the ring, in order
env = make_env(ctx, {s2: gen(ctx, 1), s12: gen(ctx, 2),
                     d: poly(ctx, {(0, 0, 0): 4, (1, 0, 0): -2})})   # d -> 4 - 2*eps
fr = walk((d - 2) / (s2 * s12), env, {}, ctx)           # exact numerator/denominator pair
num_terms, den_terms = canonical_terms(fr)

fr is a Frac holding the reduced numerator $2-2\varepsilon$ and denominator $s_2 s_{12}$ (the self-test checks exactly this identity), and canonical_terms returns each as a list of [exponents, "coefficient"] entries with exponents in the order given to mk_ctx, or None when the whole expression is exactly zero. The third argument of walk is a cache dictionary, one per expression tree, so shared subexpressions are evaluated once.

Routines

Loading a system

Transport and local solutions

Expansion in ε

ε → 0 limits from an existing grid (from wayfinder.epslimit import ..., or from wayfinder import epslimit; formerly the separate eps-extrapolator tool)

Boundary values (AMFlow conventions)

Checks, records, quadrature

Conic third-kind reduction (from wayfinder import conic_thirdkind)

Exact rational-function algebra and export (from wayfinder.flintexport import ..., or from wayfinder import flintexport; needs python-flint and SymPy; formerly the separate flint-export tool)

Connection sampler (sampler.py)

Used on this site

Requirements and source

Python 3 with mpmath. python-flint is optional (faster exact arithmetic, backend="acb" and attach_acb_fast) except for flintexport, which requires it; numpy is only a fallback root finder. sympy is needed by frob_scalar, boundary_branches, conic_thirdkind, flintexport and the sampler, and is imported only when those are used. epslimit needs nothing beyond mpmath (its gmpy2 backend helps but is optional). Importing wayfinder loads neither epslimit nor flintexport; from wayfinder import epslimit, flintexport brings them in on first use. The sampler also needs Gatekeeper, the canonical_form module of Counterweight and (for graph-defined targets) the Landau Alphabet package importable under tools/, plus an amflow_cli binary built from the amflow-cpp-dev fork (see AMFlow; AMFLOW_CLI names it, AMFLOW_IBP_CACHE sets the reduction cache). Set OPENBLAS_NUM_THREADS to 2 or less for large transports.

Tests: python3 -m pytest tests/ -q from the package directory runs the unit suite (ten to fifteen minutes), or python3 -m pytest tools/wayfinder -x -q from the repository root; the two control scripts are shown above. python3 tests/test_epslimit_synthetic.py and python3 tests/test_epslimit_fixtures.py check epslimit against expansions known in closed form in a few seconds, and python3 tests/test_flintexport.py runs the flintexport identities in under a second (skipped without python-flint); all three are also collected by pytest. python tests/test_sampler.py runs the sampler's synthetic round-trips and a negative control (one deliberately perturbed derivative value that the fit must flag) in a few seconds; WAYFINDER_RUN_AMFLOW=1 adds the AMFlow-backed box tests (ten to twenty minutes). Tests that compare against large reference systems not distributed with the package are skipped unless WAYFINDER_REFERENCE_FIXTURES points at a local copy. The code is in tools/wayfinder/ in BootLoops' bootloops-dev repository (GitHub organization BootLoops-ai), released under the MIT license.

← back to the tools index