API reference

The trusted core (rational, poly, model, certificate, checker) uses only the standard library and exact arithmetic. The search side (solver, candidate, search) uses numpy and only proposes certificates.

Checker

twindex.checker

The certificate checker (trusted core). Exact rational arithmetic only.

Mathematical basis (proved in docs/THEORY.md):

(T1) For theta >= 0 the problem has a unique classical solution u_theta. (T2) If p(0) = p(1) = 0 and |-p'' + q p^3 - 1| <= rho on [0, 1], then |u_theta(x) - p(x)| <= rho * x (1 - x) / 2 on [0, 1]. (T3) For fixed p the residual is affine in theta, so for theta in c +- h, sup |r(.; theta)| <= sup |r(.; c)| + sum_j h_j sup |B_j p^3|.

Consequences at a sensor x_i with s_i = rho * x_i (1 - x_i) / 2: |p(x_i) - y_i| + s_i <= eps proves |u(x_i) - y_i| <= eps, |p(x_i) - y_i| - s_i > eps proves |u(x_i) - y_i| > eps.

check_witness(problem, theta, p, depth=DEFAULT_DEPTH)

Claim: theta is in C_y. VALID proves it; INVALID means malformed or disproved.

Source code in twindex/checker.py

def check_witness(
    problem: Problem, theta: Sequence[Fraction | int | str], p: Poly, depth: int = DEFAULT_DEPTH
) -> Report:
    """Claim: theta is in C_y. VALID proves it; INVALID means malformed or disproved."""
    bad = _bad_inputs(problem, p)
    if bad:
        return _invalid(bad)
    try:
        t = tuple(parse_rational(v) for v in theta)
    except (TypeError, ValueError) as e:
        return _invalid(f"theta: {e}")
    if len(t) != problem.num_params:
        return _invalid(f"theta must have {problem.num_params} entries")
    if not problem.contains(t):
        return _invalid("theta lies outside the prior box", theta=[_fmt(v) for v in t])
    bc = _boundary_ok(p, problem.model.boundary)
    if bc:
        return _invalid(bc)
    rho = sup_abs_bound(problem.model.residual(p, t), depth)
    all_inside, outside, rows = _sensor_analysis(problem, p, rho)
    details = {"theta": [_fmt(v) for v in t], "rho": _fmt(rho), "sensors": rows}
    if outside:
        return Report(
            Verdict.INVALID,
            tuple(f"proved incompatible at sensor x = {rows[i]['x']}" for i in outside),
            details,
            code="INVALID_FALSE_CLAIM",
        )
    if all_inside:
        return Report(Verdict.VALID, ("theta is data-compatible",), details)
    return Report(
        Verdict.INCONCLUSIVE,
        ("residual bound too large to decide every sensor; refine p or raise depth",),
        details,
    )

check_pair(problem, w1, w2, claimed_gap, depth=DEFAULT_DEPTH)

Claim: both thetas are in C_y and |theta_k - theta'_k| >= claimed_gap.

Source code in twindex/checker.py

def check_pair(
    problem: Problem,
    w1: tuple[Sequence[Fraction | int | str], Poly],
    w2: tuple[Sequence[Fraction | int | str], Poly],
    claimed_gap: Fraction | int | str,
    depth: int = DEFAULT_DEPTH,
) -> Report:
    """Claim: both thetas are in C_y and |theta_k - theta'_k| >= claimed_gap."""
    try:
        gap_claim = parse_rational(claimed_gap)
        t1 = tuple(parse_rational(v) for v in w1[0])
        t2 = tuple(parse_rational(v) for v in w2[0])
    except (TypeError, ValueError) as e:
        return _invalid(f"pair: {e}")
    if gap_claim < 0:
        return _invalid("claimed gap must be nonnegative")
    r1 = check_witness(problem, *w1, depth=depth)
    r2 = check_witness(problem, *w2, depth=depth)
    if not isinstance(problem, Problem):
        return r1
    k = problem.target
    details: dict[str, Any] = {"witnesses": [r1.to_dict(), r2.to_dict()]}
    if len(t1) == problem.num_params and len(t2) == problem.num_params:
        actual = abs(t1[k] - t2[k])
        details["target_index"] = k
        details["actual_gap"] = _fmt(actual)
        details["claimed_gap"] = _fmt(gap_claim)
        if actual < gap_claim:
            return Report(
                Verdict.INVALID,
                (f"claimed gap {_fmt(gap_claim)} exceeds the exact gap {_fmt(actual)}",),
                details,
                code="INVALID_FALSE_CLAIM",
            )
    verdicts = {r1.verdict, r2.verdict}
    for worst in (Verdict.INVALID, Verdict.INCONCLUSIVE):
        if worst in verdicts:
            reasons = tuple(
                f"witness {i + 1}: {msg}"
                for i, r in enumerate((r1, r2))
                if r.verdict is worst
                for msg in r.reasons
            )
            return Report(worst, reasons, details)
    diameter = abs(t1[k] - t2[k])
    details["diameter_lower_bound"] = _fmt(diameter)
    details["recovery_error_lower_bound"] = _fmt(diameter / 2)
    return Report(
        Verdict.VALID,
        (
            f"both witnesses are data-compatible; D_Q(y) >= {_fmt(diameter)}, so every "
            f"estimator of theta_{k} has worst-case error >= {_fmt(diameter / 2)} "
            "on these data",
        ),
        details,
    )

box_residual_bound(center, half_width, p, depth=DEFAULT_DEPTH, cache=None, model=V1_MODEL)

Upper bound of sup_{x in [0,1], theta in center +- half_width} |r(x; theta)| (T3).

Source code in twindex/checker.py

def box_residual_bound(
    center: Sequence[Fraction],
    half_width: Sequence[Fraction],
    p: Poly,
    depth: int = DEFAULT_DEPTH,
    cache: _CandidateCache | None = None,
    model: Model = V1_MODEL,
) -> Fraction:
    """Upper bound of sup_{x in [0,1], theta in center +- half_width} |r(x; theta)| (T3)."""
    cc = cache if cache is not None else _CandidateCache(p, depth, model)
    rho = sup_abs_bound(cc.residual(center), depth)
    for j, h in enumerate(half_width):
        if h:
            rho += h * cc.sup_bj_fp(j)
    return rho

monotone_hypotheses(problem, depth=12)

Check the hypotheses of THEORY Theorem 6: alpha, beta >= 0 and g >= 0 on [0, 1].

VALID: proved. INVALID: disproved. INCONCLUSIVE: g >= 0 not decided at depth.

Source code in twindex/checker.py

def monotone_hypotheses(problem: Problem, depth: int = 12) -> tuple[Verdict, str]:
    """Check the hypotheses of THEORY Theorem 6: alpha, beta >= 0 and g >= 0 on [0, 1].

    VALID: proved. INVALID: disproved. INCONCLUSIVE: g >= 0 not decided at ``depth``.
    """
    alpha, beta = problem.model.boundary
    if alpha < 0 or beta < 0:
        return Verdict.INVALID, "monotone boxes need boundary values >= 0"
    status = nonneg_status(problem.model.source, depth)
    if status is False:
        return Verdict.INVALID, "monotone boxes need a source g >= 0; g is negative somewhere"
    if status is None:
        return Verdict.INCONCLUSIVE, "could not certify g >= 0 for monotone boxes"
    return Verdict.VALID, "monotone hypotheses hold"

check_box(problem, center, half_width, p, status, depth=DEFAULT_DEPTH, sensitivities=(), corners=None)

Claim: B = center +- half_width lies inside C_y ('inside') or misses it ('outside').

With sensitivities the candidate is the affine family p + sum_j D_j s_j. With corners = (p_lo, p_hi) (and p=None) the box is decided by the monotone two-corner bound of THEORY Theorem 6.

Source code in twindex/checker.py

def check_box(
    problem: Problem,
    center: Sequence[Fraction | int | str],
    half_width: Sequence[Fraction | int | str],
    p: Poly | None,
    status: str,
    depth: int = DEFAULT_DEPTH,
    sensitivities: Sequence[Poly] = (),
    corners: tuple[Poly, Poly] | None = None,
) -> Report:
    """Claim: B = center +- half_width lies inside C_y ('inside') or misses it ('outside').

    With ``sensitivities`` the candidate is the affine family p + sum_j D_j s_j.
    With ``corners = (p_lo, p_hi)`` (and ``p=None``) the box is decided by the
    monotone two-corner bound of THEORY Theorem 6.
    """
    if corners is not None:
        return _check_monotone_box(problem, center, half_width, p, status, depth,
                                   sensitivities, corners)
    bad = _bad_inputs(problem, p)
    if bad:
        return _invalid(bad)
    assert p is not None
    try:
        c = tuple(parse_rational(v) for v in center)
        h = tuple(parse_rational(v) for v in half_width)
    except (TypeError, ValueError) as e:
        return _invalid(f"box: {e}")
    n = problem.num_params
    if len(c) != n or len(h) != n:
        return _invalid(f"center and half_width must have {n} entries")
    if any(v < 0 for v in h):
        return _invalid("half_width entries must be nonnegative")
    if status not in ("inside", "outside"):
        return _invalid("status must be 'inside' or 'outside'")
    for j, ((lo, hi), cj, hj) in enumerate(zip(problem.box, c, h, strict=True)):
        if cj - hj < lo or cj + hj > hi:
            return _invalid(f"box exceeds the prior box in coordinate {j}")
    bc = _boundary_ok(p, problem.model.boundary)
    if bc:
        return _invalid(bc)
    sens = tuple(sensitivities)
    bad_sens = _check_sensitivities(sens, n)
    if bad_sens:
        return _invalid(bad_sens)
    try:
        rho = _family_bound(problem.model, c, h, p, sens, depth, None)
    except ExpansionTooLarge as e:
        return _invalid(str(e))
    all_inside, outside, rows = _sensor_analysis(problem, p, rho, sens, h)
    details = {
        "center": [_fmt(v) for v in c],
        "half_width": [_fmt(v) for v in h],
        "rho": _fmt(rho),
        "sensors": rows,
    }
    all_msg = "every theta in the box is data-compatible"
    none_msg = "no theta in the box is data-compatible"
    undecided = Report(
        Verdict.INCONCLUSIVE, ("could not decide the box; split it or refine p",), details
    )
    if status == "inside":
        if all_inside:
            return Report(Verdict.VALID, (all_msg,), details)
        if outside:
            return Report(Verdict.INVALID, (f"proved: {none_msg}",), details,
                          code="INVALID_FALSE_CLAIM")
        return undecided
    if outside:
        return Report(Verdict.VALID, (none_msg,), details)
    if all_inside:
        return Report(Verdict.INVALID, (f"proved: {all_msg}",), details, code="INVALID_FALSE_CLAIM")
    return undecided

check_paving(problem, paving, depth=DEFAULT_DEPTH)

Verify a SIVIA paving and its claim lower <= D_Q(y) <= upper (DESIGN.md section 6).

Source code in twindex/checker.py

def check_paving(problem: Problem, paving: Paving, depth: int = DEFAULT_DEPTH) -> Report:
    """Verify a SIVIA paving and its claim lower <= D_Q(y) <= upper (DESIGN.md section 6)."""
    if not isinstance(problem, Problem):
        return _invalid("problem must be a twindex.model.Problem")
    if not isinstance(paving, Paving):
        return _invalid("paving must be a twindex.certificate.Paving")
    if paving.transform is not None:
        return _check_oriented_paving(problem, paving, depth)
    k = problem.target
    caches: dict[int, _CandidateCache] = {}
    monotone_ok: tuple[Verdict, str] | None = None
    corner_rho: dict[tuple[int, tuple[Fraction, ...]], Fraction] = {}
    counts = {"inside": 0, "outside": 0, "unresolved": 0}
    inside_lo: Fraction | None = None
    inside_hi: Fraction | None = None
    outer_lo: Fraction | None = None
    outer_hi: Fraction | None = None
    undecided: list[str] = []
    stack: list[tuple[PavingLeaf | PavingSplit, tuple[tuple[Fraction, Fraction], ...], str]] = [
        (paving.root, problem.box, "paving")
    ]
    while stack:
        node, box, where = stack.pop()
        if isinstance(node, PavingSplit):
            if node.axis >= problem.num_params:
                return _invalid(f"{where}: split axis {node.axis} exceeds the model's parameters")
            lo, hi = box[node.axis]
            if not (lo < node.at < hi):
                return _invalid(f"{where}: split point {_fmt(node.at)} not strictly inside "
                                f"[{_fmt(lo)}, {_fmt(hi)}] on axis {node.axis}")
            left = tuple((lo_, node.at) if j == node.axis else (lo_, hi_)
                         for j, (lo_, hi_) in enumerate(box))
            right = tuple((node.at, hi_) if j == node.axis else (lo_, hi_)
                          for j, (lo_, hi_) in enumerate(box))
            stack.append((node.right, right, f"{where}.children[1]"))
            stack.append((node.left, left, f"{where}.children[0]"))
            continue
        counts[node.status] += 1
        lo_k, hi_k = box[k]
        if node.status in ("inside", "unresolved"):
            outer_lo = lo_k if outer_lo is None else min(outer_lo, lo_k)
            outer_hi = hi_k if outer_hi is None else max(outer_hi, hi_k)
        if node.status == "inside":
            inside_lo = lo_k if inside_lo is None else min(inside_lo, lo_k)
            inside_hi = hi_k if inside_hi is None else max(inside_hi, hi_k)
        if node.status == "unresolved":
            continue
        if node.p is None and node.corners is None:
            return _invalid(f"{where}: a candidate-free leaf is only valid in an oriented "
                            "paving whose leaf misses the prior")
        if node.corners is not None:
            if monotone_ok is None:
                monotone_ok = monotone_hypotheses(problem)
            if monotone_ok[0] is not Verdict.VALID:
                if monotone_ok[0] is Verdict.INVALID:
                    return _invalid(f"{where}: {monotone_ok[1]}")
                undecided.append(where)
                continue
            lo_c = tuple(a for a, _ in box)
            hi_c = tuple(b for _, b in box)
            ends: list[tuple[Fraction, Fraction]] = []
            for idx, theta, sign in ((node.corners[0], lo_c, 1), (node.corners[1], hi_c, -1)):
                cand = paving.candidates[idx]
                bc = _boundary_ok(cand, problem.model.boundary)
                if bc:
                    return _invalid(f"{where}: candidates[{idx}]: {bc}")
                key = (idx, theta)
                rho = corner_rho.get(key)
                if rho is None:
                    rho = corner_rho[key] = sup_abs_bound(
                        problem.model.residual(cand, theta), depth
                    )
                ends.append((rho, Fraction(sign)))
            (rho_lo, _), (rho_hi, _) = ends
            p_lo = paving.candidates[node.corners[0]]
            p_hi = paving.candidates[node.corners[1]]
            m_inside, m_outside = True, False
            for x, y in zip(problem.sensors, problem.data, strict=True):
                low = p_hi(x) - rho_hi * barrier(x)
                high = p_lo(x) + rho_lo * barrier(x)
                if low < y - problem.epsilon or high > y + problem.epsilon:
                    m_inside = False
                if low > y + problem.epsilon or high < y - problem.epsilon:
                    m_outside = True
            if node.status == "inside":
                if m_outside:
                    return _invalid(f"{where}: labelled inside but proved disjoint from C_y",
                                    code="INVALID_FALSE_CLAIM")
                if not m_inside:
                    undecided.append(where)
            else:
                if m_inside:
                    return _invalid(f"{where}: labelled outside but proved inside C_y",
                                    code="INVALID_FALSE_CLAIM")
                if not m_outside:
                    undecided.append(where)
            continue
        assert node.p is not None
        bc = _boundary_ok(paving.candidates[node.p], problem.model.boundary)
        if bc:
            return _invalid(f"{where}: candidates[{node.p}]: {bc}")
        cache = caches.get(node.p)
        if cache is None:
            cache = caches[node.p] = _CandidateCache(
                paving.candidates[node.p], depth, problem.model
            )
        center = tuple((lo + hi) / 2 for lo, hi in box)
        half = tuple((hi - lo) / 2 for lo, hi in box)
        sens = tuple(paving.candidates[j] for j in (node.s or ()))
        bad_sens = _check_sensitivities(sens, problem.num_params)
        if bad_sens:
            return _invalid(f"{where}: {bad_sens}")
        try:
            rho = _family_bound(problem.model, center, half, cache.p, sens, depth, cache)
        except ExpansionTooLarge as e:
            return _invalid(f"{where}: {e}")
        all_inside, outside, _ = _sensor_analysis(problem, cache.p, rho, sens, half)
        if node.status == "inside":
            if outside:
                return _invalid(f"{where}: labelled inside but proved disjoint from C_y",
                                code="INVALID_FALSE_CLAIM")
            if not all_inside:
                undecided.append(where)
        else:
            if all_inside:
                return _invalid(f"{where}: labelled outside but proved inside C_y",
                                code="INVALID_FALSE_CLAIM")
            if not outside:
                undecided.append(where)

    lower_cert = None
    if inside_lo is not None and inside_hi is not None:
        lower_cert = inside_hi - inside_lo
    upper_cert = Fraction(0)
    if outer_lo is not None and outer_hi is not None:
        upper_cert = outer_hi - outer_lo
    details: dict[str, Any] = {
        "leaves": counts,
        "target_index": k,
        "certified_lower": None if lower_cert is None else _fmt(lower_cert),
        "certified_upper": _fmt(upper_cert),
        "claimed_lower": None if paving.lower is None else _fmt(paving.lower),
        "claimed_upper": _fmt(paving.upper),
        "compatible_set_empty": outer_lo is None,
    }
    if paving.lower is not None:
        if lower_cert is None:
            return Report(Verdict.INVALID, ("claim.lower needs at least one inside leaf",), details)
        if paving.lower > lower_cert:
            return Report(
                Verdict.INVALID,
                (f"claim.lower {_fmt(paving.lower)} exceeds the certified {_fmt(lower_cert)}",),
                details,
                code="INVALID_FALSE_CLAIM",
            )
    if paving.upper < upper_cert:
        return Report(
            Verdict.INVALID,
            (f"claim.upper {_fmt(paving.upper)} is below the certified {_fmt(upper_cert)}",),
            details,
            code="INVALID_FALSE_CLAIM",
        )
    if undecided:
        details["undecided_leaves"] = undecided[:20]
        return Report(
            Verdict.INCONCLUSIVE,
            (f"{len(undecided)} labelled leaf/leaves could not be verified",),
            details,
        )
    if outer_lo is None:
        msg = "every leaf is outside: C_y is empty"
    else:
        lo_txt = "-" if lower_cert is None else _fmt(lower_cert)
        msg = f"certified {lo_txt} <= D_Q(y) <= {_fmt(upper_cert)} for theta_{k}"
    return Report(Verdict.VALID, (msg,), details)

check_certificate(cert, depth=DEFAULT_DEPTH)

Check a parsed Certificate or a raw decoded-JSON dict.

Source code in twindex/checker.py

def check_certificate(cert: Certificate | dict[str, Any], depth: int = DEFAULT_DEPTH) -> Report:
    """Check a parsed Certificate or a raw decoded-JSON dict."""
    if isinstance(cert, dict) and cert.get("kind") in ("box_ode", "paving_ode"):
        from .ode_bounds import check_certificate_bounds_ode

        return check_certificate_bounds_ode(cert, depth)
    if isinstance(cert, dict) and cert.get("kind") in ("witness_ode", "pair_ode"):
        from .ode import check_certificate_ode  # polynomial ODE systems (THEORY 7.3)

        return check_certificate_ode(cert, depth)
    if isinstance(cert, dict) and cert.get("kind") in ("witness_heat", "pair_heat"):
        from .parabolic import check_certificate_heat

        return check_certificate_heat(cert, min(depth, 4))
    if isinstance(cert, dict) and cert.get("kind") in ("witness2", "pair2"):
        from .elliptic2d import check_certificate2  # 2-D certificates (THEORY 7.1)

        return check_certificate2(cert, min(depth, 4))
    if not isinstance(cert, Certificate):
        try:
            cert = parse_certificate(cert)
        except CertificateError as e:
            unsupported = (isinstance(cert, dict) and isinstance(cert.get("kind"), str)
                           and cert["kind"] not in {"witness", "pair", "box", "paving"})
            return _invalid(f"malformed certificate: {e}",
                            code="UNSUPPORTED_MODEL" if unsupported else "INVALID_MALFORMED")
    if cert.kind == "witness":
        w = cert.witnesses[0]
        return check_witness(cert.problem, w.theta, w.p, depth)
    if cert.kind == "pair":
        a, b = cert.witnesses
        assert cert.gap is not None
        return check_pair(cert.problem, (a.theta, a.p), (b.theta, b.p), cert.gap, depth)
    if cert.kind == "paving":
        assert cert.paving is not None
        return check_paving(cert.problem, cert.paving, depth)
    assert cert.box is not None
    bx = cert.box
    return check_box(
        cert.problem, bx.center, bx.half_width, bx.p, bx.status, depth, bx.sensitivities,
        bx.corners,
    )

Certificates

twindex.certificate

Certificate JSON format v1: strict parsing, validation and construction (trusted core).

See docs/DESIGN.md section 6. Every number is a JSON string holding an exact rational. Unknown keys are rejected so that typos cannot silently drop a field.

CertificateError

Bases: ValueError

The certificate is malformed (a structural or hypothesis violation).

Source code in twindex/certificate.py

class CertificateError(ValueError):
    """The certificate is malformed (a structural or hypothesis violation)."""

Transform dataclass

Oriented paving coordinates: theta = origin + matrix @ phi, phi in the root box.

Source code in twindex/certificate.py

@dataclass(frozen=True)
class Transform:
    """Oriented paving coordinates: theta = origin + matrix @ phi, phi in the root box."""

    origin: tuple[Fraction, ...]
    matrix: tuple[tuple[Fraction, ...], ...]
    root: tuple[tuple[Fraction, Fraction], ...]
    inverse: tuple[tuple[Fraction, ...], ...] = ()

    def __post_init__(self) -> None:
        origin = _typed_rationals(self.origin, "transform.origin")
        n = len(origin)
        if n == 0:
            raise CertificateError("transform.origin must be nonempty")
        try:
            matrix = tuple(_typed_rationals(row, f"transform.matrix[{i}]", n)
                           for i, row in enumerate(self.matrix))
        except TypeError:
            raise CertificateError("transform.matrix must be a list of rows") from None
        if len(matrix) != n:
            raise CertificateError(f"transform.matrix must be {n} x {n}")
        root = []
        for i, iv in enumerate(self.root):
            lo, hi = _typed_rationals(iv, f"transform.root[{i}]", 2)
            if lo > hi:
                raise CertificateError("transform.root intervals must satisfy lo <= hi")
            root.append((lo, hi))
        if len(root) != n:
            raise CertificateError(f"transform.root must have {n} intervals")
        try:
            inverse = exact_inverse(matrix)
        except SingularMatrixError as e:
            raise CertificateError(f"transform.matrix: {e}") from None
        object.__setattr__(self, "origin", origin)
        object.__setattr__(self, "matrix", matrix)
        object.__setattr__(self, "root", tuple(root))
        object.__setattr__(self, "inverse", inverse)

parse_certificate(obj)

Validate a decoded JSON object and return a typed Certificate.

Raises CertificateError on any structural problem.

Source code in twindex/certificate.py

def parse_certificate(obj: Any) -> Certificate:
    """Validate a decoded JSON object and return a typed Certificate.

    Raises CertificateError on any structural problem.
    """
    m = _require_mapping(obj, "certificate")
    if m.get("format") != FORMAT:
        raise CertificateError(f"format must be {FORMAT!r}")
    version = m.get("version")
    if type(version) is not int or version not in SUPPORTED_VERSIONS:
        raise CertificateError(
            f"unsupported version {version!r}; this checker reads {SUPPORTED_VERSIONS}"
        )
    kind = m.get("kind")
    if kind not in KINDS:
        raise CertificateError(f"kind must be one of {KINDS}, got {kind!r}")
    required = {"format", "version", "kind", "problem"}
    if kind == "witness":
        required |= {"witnesses"}
        allowed = required | {"provenance"}
    elif kind == "pair":
        required |= {"witnesses", "claim"}
        allowed = required | {"provenance"}
    elif kind == "box":
        required |= {"box", "claim"}
        allowed = required | {"provenance"}
    else:  # paving
        required |= {"candidates", "paving", "claim"}
        allowed = required | {"provenance", "transform"}
    _check_keys(m, allowed, required, "certificate")
    problem = _problem(m["problem"], version)
    provenance = m.get("provenance", {})
    if not isinstance(provenance, Mapping):
        raise CertificateError("provenance must be a JSON object")

    if kind in ("witness", "pair"):
        ws = m["witnesses"]
        expected = 1 if kind == "witness" else 2
        if not isinstance(ws, list) or len(ws) != expected:
            raise CertificateError(f"kind {kind!r} needs exactly {expected} witness(es)")
        witnesses = tuple(_witness(w, f"witnesses[{i}]") for i, w in enumerate(ws))
        gap = None
        if kind == "pair":
            claim = _require_mapping(m["claim"], "claim")
            _check_keys(claim, {"gap"}, {"gap"}, "claim")
            gap = _rational(claim["gap"], "claim.gap")
            if gap < 0:
                raise CertificateError("claim.gap must be nonnegative")
        return Certificate(kind, problem, witnesses, gap, None, provenance)

    if kind == "paving":
        paving = _paving(m["candidates"], m["paving"], m["claim"])
        if "transform" in m:
            paving = Paving(paving.candidates, paving.root, paving.lower, paving.upper,
                            _transform(m["transform"]))
        return Certificate(kind, problem, (), None, None, provenance, paving)

    b = _require_mapping(m["box"], "box")
    _check_keys(b, _BOX_KEYS, {"center", "half_width"}, "box")
    if ("p" in b) == ("corners" in b):
        raise CertificateError("box needs exactly one of 'p' and 'corners'")
    corners = None
    if "corners" in b:
        if "sensitivities" in b:
            raise CertificateError("a corners box takes no sensitivities")
        cs = b["corners"]
        if not isinstance(cs, list) or len(cs) != 2:
            raise CertificateError("box.corners must be [p_lower_corner, p_upper_corner]")
        corners = (_poly(cs[0], "box.corners[0]"), _poly(cs[1], "box.corners[1]"))
    claim = _require_mapping(m["claim"], "claim")
    _check_keys(claim, {"status"}, {"status"}, "claim")
    status = claim["status"]
    if status not in ("inside", "outside"):
        raise CertificateError("claim.status must be 'inside' or 'outside'")
    box = BoxClaim(
        center=_rational_list(b["center"], "box.center"),
        half_width=_rational_list(b["half_width"], "box.half_width"),
        p=_poly(b["p"], "box.p") if "p" in b else None,
        status=status,
        sensitivities=tuple(
            _poly(s, f"box.sensitivities[{i}]") for i, s in enumerate(_list(b, "sensitivities"))
        ),
        corners=corners,
    )
    if any(h < 0 for h in box.half_width):
        raise CertificateError("box.half_width entries must be nonnegative")
    return Certificate(kind, problem, (), None, box, provenance)

loads_certificate_json(text)

Decode certificate JSON text under strict budgets.

Rejected: JSON floats, NaN and Infinity; duplicate object keys (ambiguous certificates); integer literals longer than MAX_JSON_INT_DIGITS digits; nesting deeper than MAX_JSON_DEPTH; input larger than MAX_CERTIFICATE_BYTES.

Source code in twindex/certificate.py

def loads_certificate_json(text: str) -> Any:
    """Decode certificate JSON text under strict budgets.

    Rejected: JSON floats, NaN and Infinity; duplicate object keys (ambiguous
    certificates); integer literals longer than MAX_JSON_INT_DIGITS digits;
    nesting deeper than MAX_JSON_DEPTH; input larger than MAX_CERTIFICATE_BYTES.
    """
    if not isinstance(text, str):
        raise CertificateError("certificate text must be a string")
    if len(text.encode("utf-8")) > MAX_CERTIFICATE_BYTES:
        raise CertificateError(f"certificate larger than {MAX_CERTIFICATE_BYTES} bytes")
    if not _json_depth_ok(text, MAX_JSON_DEPTH):
        raise CertificateError(f"JSON nesting deeper than {MAX_JSON_DEPTH}")

    def no_floats(s: str) -> Any:
        raise CertificateError(f"JSON float {s[:40]} rejected: write rationals as strings")

    def bounded_int(s: str) -> int:
        if len(s.lstrip("-")) > MAX_JSON_INT_DIGITS:
            raise CertificateError("JSON integer literal too long")
        return int(s)

    def unique_keys(pairs: list[tuple[str, Any]]) -> dict[str, Any]:
        out: dict[str, Any] = {}
        for k, v in pairs:
            if k in out:
                raise CertificateError(f"duplicate JSON key {k[:40]!r}")
            out[k] = v
        return out

    try:
        return json.loads(text, parse_float=no_floats, parse_constant=no_floats,
                          parse_int=bounded_int, object_pairs_hook=unique_keys)
    except json.JSONDecodeError as e:
        raise CertificateError(f"invalid JSON: {e}") from None
    except RecursionError:
        raise CertificateError("JSON nesting too deep") from None

read_certificate_json(path)

Decode a certificate file under the budgets of :func:loads_certificate_json.

Source code in twindex/certificate.py

def read_certificate_json(path: str | Path) -> Any:
    """Decode a certificate file under the budgets of :func:`loads_certificate_json`."""
    p = Path(path)
    if p.stat().st_size > MAX_CERTIFICATE_BYTES:
        raise CertificateError(f"certificate larger than {MAX_CERTIFICATE_BYTES} bytes")
    with open(p, encoding="utf-8") as f:
        return loads_certificate_json(f.read())

load_certificate(path)

Read and parse a (1-D) certificate file. JSON numbers that are floats are rejected.

Source code in twindex/certificate.py

def load_certificate(path: str | Path) -> Certificate:
    """Read and parse a (1-D) certificate file. JSON numbers that are floats are rejected."""
    return parse_certificate(read_certificate_json(path))

make_box_certificate(problem, center, half_width, p, status, provenance=None, sensitivities=None, corners=None)

Box certificate with a single candidate p (optionally an affine family via sensitivities) or, with p=None, a monotone corners=(p_lo, p_hi) box.

Source code in twindex/certificate.py

def make_box_certificate(
    problem: Problem,
    center: Sequence[Fraction],
    half_width: Sequence[Fraction],
    p: Poly | None,
    status: str,
    provenance: Mapping[str, Any] | None = None,
    sensitivities: Sequence[Poly] | None = None,
    corners: tuple[Poly, Poly] | None = None,
) -> dict[str, Any]:
    """Box certificate with a single candidate ``p`` (optionally an affine family via
    ``sensitivities``) or, with ``p=None``, a monotone ``corners=(p_lo, p_hi)`` box."""
    d = _envelope("box", problem, provenance)
    d["box"] = {
        "center": _fmt_list(center),
        "half_width": _fmt_list(half_width),
    }
    if p is not None:
        d["box"]["p"] = _fmt_list(p.coeffs or (Fraction(0),))
    if corners is not None:
        d["box"]["corners"] = [_fmt_list(c.coeffs or (Fraction(0),)) for c in corners]
    if sensitivities:
        d["box"]["sensitivities"] = [_fmt_list(s.coeffs or (Fraction(0),)) for s in sensitivities]
    d["claim"] = {"status": status}
    return _reorder(d)

make_paving_certificate(problem, candidates, tree, upper, lower=None, provenance=None, transform=None)

Build a paving certificate.

tree uses the JSON node shape of DESIGN.md section 6; at values may be Fractions and are formatted here. transform=(origin, matrix, root) makes it oriented (THEORY Theorem 8): leaves are boxes in phi, theta = origin + matrix @ phi.

Source code in twindex/certificate.py

def make_paving_certificate(
    problem: Problem,
    candidates: Sequence[Poly],
    tree: Mapping[str, Any],
    upper: Fraction,
    lower: Fraction | None = None,
    provenance: Mapping[str, Any] | None = None,
    transform: tuple[Sequence[Fraction], Sequence[Sequence[Fraction]],
                     Sequence[Sequence[Fraction]]] | None = None,
) -> dict[str, Any]:
    """Build a paving certificate.

    ``tree`` uses the JSON node shape of DESIGN.md section 6; ``at`` values may be
    Fractions and are formatted here. ``transform=(origin, matrix, root)`` makes it
    oriented (THEORY Theorem 8): leaves are boxes in phi, theta = origin + matrix @ phi.
    """
    d = _envelope("paving", problem, provenance)
    d["candidates"] = [_fmt_list(p.coeffs or (Fraction(0),)) for p in candidates]
    d["paving"] = _format_tree(tree)
    claim: dict[str, str] = {}
    if lower is not None:
        claim["lower"] = format_rational(lower)
    claim["upper"] = format_rational(upper)
    d["claim"] = claim
    if transform is not None:
        origin, matrix, root = transform
        d["transform"] = {
            "origin": _fmt_list(origin),
            "matrix": [_fmt_list(row) for row in matrix],
            "root": [_fmt_list(iv) for iv in root],
        }
    return _reorder(d)

dump_certificate(cert, path)

Write a certificate dict as pretty JSON (validating it first).

Source code in twindex/certificate.py

def dump_certificate(cert: Mapping[str, Any], path: str | Path) -> None:
    """Write a certificate dict as pretty JSON (validating it first)."""
    parse_certificate(cert)
    with open(path, "w", encoding="utf-8") as f:
        json.dump(cert, f, indent=2)
        f.write("\n")

Model

twindex.model

The boundary value problem and its exact building blocks (trusted core).

-u''(x) + q(x; theta) f(u(x)) = g(x) on (0, 1),   u(0) = alpha, u(1) = beta,
q(x; theta) = sum_{j=0..d} theta_j B_j(x),   B_j(x) = C(d, j) x^j (1 - x)^(d - j),
f(u) = sum_m a_m u^m   with odd m and a_m >= 0,   g a polynomial.

The default Model() is the original v1 problem: d = 3, f(u) = u^3, g = 1. Every hypothesis the checker relies on (theta >= 0 so that q >= 0, odd powers with nonnegative coefficients so that f is nondecreasing) is enforced here.

ProblemError

Bases: ValueError

A problem instance or parameter vector violates a hypothesis.

Source code in twindex/model.py

class ProblemError(ValueError):
    """A problem instance or parameter vector violates a hypothesis."""

Model dataclass

The equation: coefficient degree, source term and nonlinearity.

nonlinearity lists (m, a_m) pairs; zero coefficients are dropped and the pairs are stored sorted by power.

Source code in twindex/model.py

@dataclass(frozen=True)
class Model:
    """The equation: coefficient degree, source term and nonlinearity.

    ``nonlinearity`` lists ``(m, a_m)`` pairs; zero coefficients are dropped and
    the pairs are stored sorted by power.
    """

    degree: int = 3
    source: Poly = field(default_factory=lambda: Poly([1]))
    nonlinearity: tuple[tuple[int, Fraction], ...] = ((3, Fraction(1)),)
    boundary: tuple[Fraction, Fraction] = (Fraction(0), Fraction(0))  # (u(0), u(1))

    def __post_init__(self) -> None:
        if isinstance(self.degree, bool) or not isinstance(self.degree, int):
            raise ProblemError("coefficient degree must be an int")
        if not 0 <= self.degree <= MAX_COEFFICIENT_DEGREE:
            raise ProblemError(f"coefficient degree must be in 0..{MAX_COEFFICIENT_DEGREE}")
        if not isinstance(self.source, Poly):
            raise ProblemError("source must be a Poly")
        if self.source.degree() > MAX_SOURCE_DEGREE:
            raise ProblemError(f"source degree must be at most {MAX_SOURCE_DEGREE}")
        object.__setattr__(self, "nonlinearity", _nonlinearity(self.nonlinearity))
        try:
            left, right = self.boundary
            bc = (parse_rational(left), parse_rational(right))
        except (TypeError, ValueError) as e:
            raise ProblemError(f"boundary must be two exact rationals: {e}") from None
        object.__setattr__(self, "boundary", bc)

    @property
    def num_params(self) -> int:
        return self.degree + 1

    @cached_property
    def basis(self) -> tuple[Poly, ...]:
        return bernstein_basis(self.degree)

    def theta(self, theta: Sequence[Fraction | int | str]) -> tuple[Fraction, ...]:
        """Parse a parameter vector exactly and check its length."""
        try:
            t = tuple(parse_rational(v) for v in theta)
        except (TypeError, ValueError) as e:
            raise ProblemError(f"theta: {e}") from None
        if len(t) != self.num_params:
            raise ProblemError(f"theta must have {self.num_params} entries, got {len(t)}")
        return t

    def q(self, theta: Sequence[Fraction | int | str]) -> Poly:
        """q(x; theta) = sum_j theta_j B_j(x)."""
        q = Poly()
        for tj, bj in zip(self.theta(theta), self.basis, strict=True):
            q = q + bj.scale(tj)
        return q

    def f_of(self, p: Poly) -> Poly:
        """f(p) = sum_m a_m p^m, exactly."""
        out = Poly()
        for m, a in self.nonlinearity:
            out = out + (p**m).scale(a)
        return out

    def residual(self, p: Poly, theta: Sequence[Fraction | int | str]) -> Poly:
        """r = -p'' + q(.; theta) f(p) - g, exactly."""
        return -p.deriv(2) + self.q(theta) * self.f_of(p) - self.source

    @property
    def is_v1(self) -> bool:
        return self == V1_MODEL

theta(theta)

Parse a parameter vector exactly and check its length.

Source code in twindex/model.py

def theta(self, theta: Sequence[Fraction | int | str]) -> tuple[Fraction, ...]:
    """Parse a parameter vector exactly and check its length."""
    try:
        t = tuple(parse_rational(v) for v in theta)
    except (TypeError, ValueError) as e:
        raise ProblemError(f"theta: {e}") from None
    if len(t) != self.num_params:
        raise ProblemError(f"theta must have {self.num_params} entries, got {len(t)}")
    return t

q(theta)

q(x; theta) = sum_j theta_j B_j(x).

Source code in twindex/model.py

def q(self, theta: Sequence[Fraction | int | str]) -> Poly:
    """q(x; theta) = sum_j theta_j B_j(x)."""
    q = Poly()
    for tj, bj in zip(self.theta(theta), self.basis, strict=True):
        q = q + bj.scale(tj)
    return q

f_of(p)

f(p) = sum_m a_m p^m, exactly.

Source code in twindex/model.py

def f_of(self, p: Poly) -> Poly:
    """f(p) = sum_m a_m p^m, exactly."""
    out = Poly()
    for m, a in self.nonlinearity:
        out = out + (p**m).scale(a)
    return out

residual(p, theta)

r = -p'' + q(.; theta) f(p) - g, exactly.

Source code in twindex/model.py

def residual(self, p: Poly, theta: Sequence[Fraction | int | str]) -> Poly:
    """r = -p'' + q(.; theta) f(p) - g, exactly."""
    return -p.deriv(2) + self.q(theta) * self.f_of(p) - self.source

Problem dataclass

A problem instance. All values are exact rationals.

Invariants (checked in __post_init__): sensors lie in the open interval (0, 1) and are distinct; one datum per sensor; epsilon > 0; the box has one interval per model parameter with 0 <= lo <= hi; 0 <= target < num_params.

Source code in twindex/model.py

@dataclass(frozen=True)
class Problem:
    """A problem instance. All values are exact rationals.

    Invariants (checked in ``__post_init__``): sensors lie in the open interval
    (0, 1) and are distinct; one datum per sensor; epsilon > 0; the box has one
    interval per model parameter with 0 <= lo <= hi; 0 <= target < num_params.
    """

    sensors: tuple[Fraction, ...]
    data: tuple[Fraction, ...]
    epsilon: Fraction
    box: tuple[tuple[Fraction, Fraction], ...]
    target: int = 0
    model: Model = V1_MODEL

    def __post_init__(self) -> None:
        # Normalise every numeric field to an exact Fraction. parse_rational rejects
        # float, bool and anything non-rational, so no float can reach the checker
        # even when a Problem is constructed directly from Python.
        if not isinstance(self.model, Model):
            raise ProblemError("model must be a twindex.model.Model")
        try:
            sensors = tuple(parse_rational(v) for v in self.sensors)
            data = tuple(parse_rational(v) for v in self.data)
            epsilon = parse_rational(self.epsilon)
            box = []
            for interval in self.box:
                lo, hi = interval
                box.append((parse_rational(lo), parse_rational(hi)))
        except (TypeError, ValueError) as e:
            raise ProblemError(f"invalid problem field: {e}") from None
        object.__setattr__(self, "sensors", sensors)
        object.__setattr__(self, "data", data)
        object.__setattr__(self, "epsilon", epsilon)
        object.__setattr__(self, "box", tuple(box))
        n = self.model.num_params
        if not self.sensors:
            raise ProblemError("at least one sensor is required")
        if len(self.data) != len(self.sensors):
            raise ProblemError("data and sensors must have the same length")
        if len(set(self.sensors)) != len(self.sensors):
            raise ProblemError("sensors must be distinct")
        for x in self.sensors:
            if not (0 < x < 1):
                raise ProblemError(f"sensor {x} is not in the open interval (0, 1)")
        if not self.epsilon > 0:
            raise ProblemError("epsilon must be positive")
        if len(self.box) != n:
            raise ProblemError(f"box must have {n} intervals (one per model parameter)")
        for lo, hi in self.box:
            if lo < 0:
                raise ProblemError("box lower bounds must be >= 0 (q >= 0 is required)")
            if lo > hi:
                raise ProblemError("box intervals must satisfy lo <= hi")
        if isinstance(self.target, bool) or not isinstance(self.target, int):
            raise ProblemError("target must be an int")
        if not (0 <= self.target < n):
            raise ProblemError(f"target must be in 0..{n - 1}")

    @property
    def num_params(self) -> int:
        return self.model.num_params

    @staticmethod
    def build(
        sensors: Sequence[Fraction | int | str],
        data: Sequence[Fraction | int | str],
        epsilon: Fraction | int | str,
        box: Sequence[Sequence[Fraction | int | str]],
        target: int = 0,
        model: Model = V1_MODEL,
    ) -> Problem:
        """Construct from loosely typed values, parsing every number exactly."""
        parsed_box = []
        for interval in box:
            if len(interval) != 2:
                raise ProblemError("each box interval must be [lo, hi]")
            parsed_box.append((parse_rational(interval[0]), parse_rational(interval[1])))
        return Problem(
            sensors=tuple(parse_rational(s) for s in sensors),
            data=tuple(parse_rational(d) for d in data),
            epsilon=parse_rational(epsilon),
            box=tuple(parsed_box),
            target=target,
            model=model,
        )

    def contains(self, theta: Sequence[Fraction | int | str]) -> bool:
        t = self.model.theta(theta)
        return all(lo <= v <= hi for v, (lo, hi) in zip(t, self.box, strict=True))

build(sensors, data, epsilon, box, target=0, model=V1_MODEL) staticmethod

Construct from loosely typed values, parsing every number exactly.

Source code in twindex/model.py

@staticmethod
def build(
    sensors: Sequence[Fraction | int | str],
    data: Sequence[Fraction | int | str],
    epsilon: Fraction | int | str,
    box: Sequence[Sequence[Fraction | int | str]],
    target: int = 0,
    model: Model = V1_MODEL,
) -> Problem:
    """Construct from loosely typed values, parsing every number exactly."""
    parsed_box = []
    for interval in box:
        if len(interval) != 2:
            raise ProblemError("each box interval must be [lo, hi]")
        parsed_box.append((parse_rational(interval[0]), parse_rational(interval[1])))
    return Problem(
        sensors=tuple(parse_rational(s) for s in sensors),
        data=tuple(parse_rational(d) for d in data),
        epsilon=parse_rational(epsilon),
        box=tuple(parsed_box),
        target=target,
        model=model,
    )

bernstein_basis(d)

Bernstein basis B_0..B_d of degree d on [0, 1]: nonnegative there, summing to 1.

Source code in twindex/model.py

def bernstein_basis(d: int) -> tuple[Poly, ...]:
    """Bernstein basis B_0..B_d of degree d on [0, 1]: nonnegative there, summing to 1."""
    return tuple((_X**j * _ONE_MINUS_X ** (d - j)).scale(comb(d, j)) for j in range(d + 1))

q_poly(theta)

q(x; theta) for the v1 model (cubic Bernstein coefficient).

Source code in twindex/model.py

def q_poly(theta: Sequence[Fraction | int | str]) -> Poly:
    """q(x; theta) for the v1 model (cubic Bernstein coefficient)."""
    return V1_MODEL.q(theta)

residual(p, theta)

r = -p'' + q p^3 - 1 for the v1 model.

Source code in twindex/model.py

def residual(p: Poly, theta: Sequence[Fraction | int | str]) -> Poly:
    """r = -p'' + q p^3 - 1 for the v1 model."""
    return V1_MODEL.residual(p, theta)

barrier(x)

x (1 - x) / 2: the barrier factor in |u - p| <= rho * x (1 - x) / 2.

Source code in twindex/model.py

def barrier(x: Fraction | int | str) -> Fraction:
    """x (1 - x) / 2: the barrier factor in |u - p| <= rho * x (1 - x) / 2."""
    xv = parse_rational(x)
    return xv * (1 - xv) / 2

Polynomials

twindex.poly

Exact rational polynomials and certified sup-norm bounds on [0, 1].

Part of the trusted core: standard library only, no floating point.

The sup bound rests on two classical facts about the Bernstein form p(x) = sum_i b_i C(n, i) x^i (1 - x)^(n - i) on [0, 1]:

1. the basis polynomials are nonnegative and sum to 1, so |p(x)| <= max_i |b_i|;
2. b_0 = p(0) and b_n = p(1), so endpoint coefficients are attained values.

De Casteljau subdivision at x = 1/2 gives the Bernstein coefficients of the two halves exactly, which tightens the bound (fact 1 applied piecewise) while fact 2 supplies a lower bound used for pruning. All arithmetic is on Python integers.

Poly

Immutable polynomial with exact rational coefficients in the monomial basis.

coeffs[k] multiplies x**k. Trailing zero coefficients are stripped, so the zero polynomial has coeffs == () and degree -1.

Source code in twindex/poly.py

class Poly:
    """Immutable polynomial with exact rational coefficients in the monomial basis.

    ``coeffs[k]`` multiplies ``x**k``. Trailing zero coefficients are stripped, so
    the zero polynomial has ``coeffs == ()`` and degree ``-1``.
    """

    __slots__ = ("_coeffs",)
    _coeffs: tuple[Fraction, ...]

    def __init__(self, coeffs: Iterable[Fraction | int | str] = ()):
        cs = [_to_fraction(c) for c in coeffs]
        while cs and cs[-1] == 0:
            cs.pop()
        object.__setattr__(self, "_coeffs", tuple(cs))

    def __setattr__(self, name: str, value: object) -> None:
        raise AttributeError("Poly is immutable")

    def __delattr__(self, name: str) -> None:
        raise AttributeError("Poly is immutable")

    # construction helpers -------------------------------------------------
    @staticmethod
    def x() -> Poly:
        return Poly([0, 1])

    @staticmethod
    def const(c: Fraction | int | str) -> Poly:
        return Poly([c])

    # basic properties -----------------------------------------------------
    @property
    def coeffs(self) -> tuple[Fraction, ...]:
        return self._coeffs

    def degree(self) -> int:
        return len(self._coeffs) - 1

    def __call__(self, x: Fraction | int | str) -> Fraction:
        xv = _to_fraction(x)
        acc = Fraction(0)
        for c in reversed(self._coeffs):
            acc = acc * xv + c
        return acc

    # arithmetic -----------------------------------------------------------
    @staticmethod
    def _coerce(other: object) -> Poly | None:
        if isinstance(other, Poly):
            return other
        if isinstance(other, (int, Fraction, str)) and not isinstance(other, bool):
            return Poly([other])
        return None

    def __add__(self, other: object) -> Poly:
        o = self._coerce(other)
        if o is None:
            return NotImplemented
        n = max(len(self._coeffs), len(o._coeffs))
        a = self._coeffs + (Fraction(0),) * (n - len(self._coeffs))
        b = o._coeffs + (Fraction(0),) * (n - len(o._coeffs))
        return Poly(x + y for x, y in zip(a, b, strict=True))

    __radd__ = __add__

    def __neg__(self) -> Poly:
        return Poly(-c for c in self._coeffs)

    def __sub__(self, other: object) -> Poly:
        o = self._coerce(other)
        if o is None:
            return NotImplemented
        return self + (-o)

    def __rsub__(self, other: object) -> Poly:
        o = self._coerce(other)
        if o is None:
            return NotImplemented
        return o - self

    def __mul__(self, other: object) -> Poly:
        o = self._coerce(other)
        if o is None:
            return NotImplemented
        if not self._coeffs or not o._coeffs:
            return Poly()
        out = [Fraction(0)] * (len(self._coeffs) + len(o._coeffs) - 1)
        for i, a in enumerate(self._coeffs):
            if a == 0:
                continue
            for j, b in enumerate(o._coeffs):
                out[i + j] += a * b
        return Poly(out)

    __rmul__ = __mul__

    def __pow__(self, n: int) -> Poly:
        if not isinstance(n, int) or isinstance(n, bool) or n < 0:
            raise ValueError("Poly exponent must be a nonnegative int")
        result = Poly([1])
        base = self
        while n:
            if n & 1:
                result = result * base
            base = base * base
            n >>= 1
        return result

    def scale(self, c: Fraction | int | str) -> Poly:
        cv = _to_fraction(c)
        return Poly(cv * a for a in self._coeffs)

    def deriv(self, k: int = 1) -> Poly:
        if k < 0:
            raise ValueError("derivative order must be nonnegative")
        p = self
        for _ in range(k):
            p = Poly(i * c for i, c in enumerate(p._coeffs) if i > 0)
        return p

    # comparison / display -------------------------------------------------
    def __eq__(self, other: object) -> bool:
        if isinstance(other, Poly):
            return self._coeffs == other._coeffs
        if isinstance(other, (int, Fraction)) and not isinstance(other, bool):
            return self._coeffs == Poly([other])._coeffs
        return NotImplemented

    def __hash__(self) -> int:
        return hash(self._coeffs)

    def __repr__(self) -> str:
        return f"Poly({[str(c) for c in self._coeffs]})"

bernstein_coeffs(p, n=None)

Exact Bernstein coefficients of p of degree n (default max(deg p, 0)) on [0, 1].

Uses b_i = sum_{k <= i} C(i, k) / C(n, k) * a_k.

Source code in twindex/poly.py

def bernstein_coeffs(p: Poly, n: int | None = None) -> list[Fraction]:
    """Exact Bernstein coefficients of ``p`` of degree ``n`` (default ``max(deg p, 0)``) on [0, 1].

    Uses ``b_i = sum_{k <= i} C(i, k) / C(n, k) * a_k``.
    """
    d = p.degree()
    if n is None:
        n = max(d, 0)
    if n < d:
        raise ValueError(f"Bernstein degree {n} is below polynomial degree {d}")
    a = list(p.coeffs) + [Fraction(0)] * (n + 1 - len(p.coeffs))
    scaled = [a[k] / comb(n, k) for k in range(n + 1)]
    return [sum((comb(i, k) * scaled[k] for k in range(i + 1)), Fraction(0)) for i in range(n + 1)]

sup_abs_bound(p, depth=6)

Certified upper bound on max_{x in [0, 1]} |p(x)|.

Exact: never smaller than the true maximum. depth bounds the number of bisection levels (at most 2**depth leaves). Pieces whose Bernstein bound cannot exceed an attained endpoint value are pruned.

Source code in twindex/poly.py

def sup_abs_bound(p: Poly, depth: int = 6) -> Fraction:
    """Certified upper bound on ``max_{x in [0, 1]} |p(x)|``.

    Exact: never smaller than the true maximum. ``depth`` bounds the number of
    bisection levels (at most ``2**depth`` leaves). Pieces whose Bernstein bound
    cannot exceed an attained endpoint value are pruned.
    """
    if depth < 0:
        raise ValueError("depth must be nonnegative")
    if p.degree() <= 0:
        return abs(p.coeffs[0]) if p.coeffs else Fraction(0)
    nums, den = _common_denominator(bernstein_coeffs(p))
    # lower bound: attained values at 0 and 1
    best_lower = max(Fraction(abs(nums[0]), den), Fraction(abs(nums[-1]), den))
    # max-heap on piece upper bound: entries (-upper, tiebreak, level, nums, den)
    counter = 0
    heap: list[tuple[Fraction, int, int, list[int], int]] = []

    def push(c: list[int], d: int, level: int) -> None:
        nonlocal counter, best_lower
        upper = Fraction(max(abs(v) for v in c), d)
        best_lower = max(best_lower, Fraction(abs(c[0]), d), Fraction(abs(c[-1]), d))
        counter += 1
        heapq.heappush(heap, (-upper, counter, level, c, d))

    push(nums, den, 0)
    final_uppers: list[Fraction] = []
    while heap:
        neg_upper, _, level, c, d = heapq.heappop(heap)
        upper = -neg_upper
        if upper <= best_lower:
            # this and every remaining piece is dominated by an attained value
            break
        if level >= depth:
            final_uppers.append(upper)
            continue
        left, right = _split_half(c)
        scale = d << (len(c) - 1)
        push(left, scale, level + 1)
        push(right, scale, level + 1)
    return max([best_lower, *final_uppers])

nonneg_status(p, depth=8)

Decide p >= 0 on [0, 1] exactly, where possible.

Returns True when proved (every Bernstein coefficient of every leaf piece is

= 0, so each piece is a convex combination of nonnegative numbers), False when disproved (some attained value b_0 or b_n of a piece is negative), and None when depth bisection levels decide neither.

Source code in twindex/poly.py

def nonneg_status(p: Poly, depth: int = 8) -> bool | None:
    """Decide ``p >= 0`` on [0, 1] exactly, where possible.

    Returns True when proved (every Bernstein coefficient of every leaf piece is
    >= 0, so each piece is a convex combination of nonnegative numbers), False
    when disproved (some attained value b_0 or b_n of a piece is negative), and
    None when ``depth`` bisection levels decide neither.
    """
    if depth < 0:
        raise ValueError("depth must be nonnegative")
    if p.degree() <= 0:
        const = p.coeffs[0] if p.coeffs else Fraction(0)
        return const >= 0
    nums, _den = _common_denominator(bernstein_coeffs(p))
    stack: list[tuple[list[int], int]] = [(nums, 0)]
    undecided = False
    while stack:
        c, level = stack.pop()
        if c[0] < 0 or c[-1] < 0:  # attained values p(left end), p(right end) of the piece
            return False
        if min(c) >= 0:
            continue
        if level >= depth:
            undecided = True
            continue
        left, right = _split_half(c)
        stack.append((left, level + 1))
        stack.append((right, level + 1))
    return None if undecided else True

Rationals

twindex.rational

Exact rational parsing and formatting for the trusted core.

Every number that enters a certificate is an exact rational written as a string: an integer ("-3"), a fraction ("22/7"), or a finite decimal ("0.125", "1e-3"). Python floats are rejected outright, because a float has usually been rounded already and silently accepting one would hide that.

RationalError

Bases: ValueError

Raised when a value is not an acceptable exact rational.

Source code in twindex/rational.py

class RationalError(ValueError):
    """Raised when a value is not an acceptable exact rational."""

parse_rational(value)

Parse an exact rational.

Accepts int, Fraction, or a string holding an integer, a/b with b > 0, or a finite decimal with optional exponent. Rejects bool, float, nan/inf and anything else.

Source code in twindex/rational.py

def parse_rational(value: str | int | Fraction) -> Fraction:
    """Parse an exact rational.

    Accepts ``int``, ``Fraction``, or a string holding an integer, ``a/b`` with
    ``b > 0``, or a finite decimal with optional exponent. Rejects ``bool``,
    ``float``, ``nan``/``inf`` and anything else.
    """
    if isinstance(value, bool):
        raise RationalError("booleans are not rationals")
    if isinstance(value, Fraction):
        return value
    if isinstance(value, int):
        return Fraction(value)
    if isinstance(value, float):
        raise RationalError(
            f"float {value!r} rejected: write exact rationals as strings such as '1/3'"
        )
    if not isinstance(value, str):
        raise RationalError(f"expected a string rational, got {type(value).__name__}")
    s = value.strip()
    if len(s) > MAX_RATIONAL_CHARS:
        raise RationalError(f"rational string longer than {MAX_RATIONAL_CHARS} characters")
    try:
        m = _FRACTION_RE.match(s)
        if m:
            den = int(m.group(2))
            if den == 0:
                raise RationalError(f"zero denominator in {value!r}")
            return Fraction(int(m.group(1)), den)
        if _DECIMAL_RE.match(s):
            exp_match = re.search(r"[eE]([+-]?\d+)$", s)
            if exp_match and abs(int(exp_match.group(1))) > MAX_EXPONENT:
                raise RationalError(f"exponent too large in {value!r}")
            return Fraction(s)  # Fraction parses finite decimals exactly
    except RationalError:
        raise
    except (ValueError, OverflowError, ZeroDivisionError) as e:
        raise RationalError(f"cannot convert {s[:40]!r}...: {e}") from None
    raise RationalError(f"not an exact rational: {value[:80]!r}")

format_rational(q)

Canonical string form: "a" for integers, "a/b" otherwise (lowest terms).

Source code in twindex/rational.py

def format_rational(q: Fraction | int) -> str:
    """Canonical string form: ``"a"`` for integers, ``"a/b"`` otherwise (lowest terms)."""
    q = Fraction(q)
    return _exact_str(q.numerator) if q.denominator == 1 else (
        f"{_exact_str(q.numerator)}/{_exact_str(q.denominator)}"
    )

Search (untrusted)

twindex.search

Untrusted, deterministic local witness search; every output is checked exactly.

This is a heuristic lower-bound search, with no claim of global optimality. Adaptive SIVIA retains every undecided box in a certified outer enclosure.

make_problem(theta_bar, sensors, eps, box, *, target=0, model=V1_MODEL, n=48, digits=12)

Round synthetic data to decimals and certify the generating parameter.

Raises RuntimeError if rounding or the candidate prevents exact certification.

Source code in twindex/search.py

def make_problem(theta_bar, sensors, eps, box, *, target=0, model=V1_MODEL, n=48, digits=12):
    """Round synthetic data to decimals and certify the generating parameter.

    Raises RuntimeError if rounding or the candidate prevents exact certification.
    """
    theta = tuple(parse_rational(v) for v in theta_bar)
    xs = tuple(parse_rational(x) for x in sensors)
    coeffs = solve_bvp(theta, n=n, model=model)
    data = _rounded(chebval(2 * np.array(xs, dtype=float) - 1, coeffs), digits)
    problem = Problem.build(xs, data, eps, box, target, model)
    report = check_witness(problem, theta, Poly(to_candidate(coeffs, boundary=model.boundary)))
    if report.verdict is not Verdict.VALID:
        raise RuntimeError(f"synthetic-data seed did not certify: {report.reasons}")
    return problem

find_pair(problem, *, seed=None, n=48, rounds=3, digits=6, depth=6)

Locally enlarge the target gap, rationalise, then certify both witnesses.

Start from a feasible seed (default: prior midpoint). Use target-coordinate steps and least-squares sensor compensation in the other coordinates. Bisection accepts only numerical sensor deviations below 0.9 epsilon. Final exact checks also require rho/8 <= epsilon/10. Failure raises; numerical feasibility alone never produces an output certificate.

Source code in twindex/search.py

def find_pair(problem, *, seed=None, n=48, rounds=3, digits=6, depth=6):
    """Locally enlarge the target gap, rationalise, then certify both witnesses.

    Start from a feasible seed (default: prior midpoint). Use target-coordinate
    steps and least-squares sensor compensation in the other coordinates.
    Bisection accepts only numerical sensor deviations below 0.9 epsilon.
    Final exact checks also require rho/8 <= epsilon/10. Failure raises;
    numerical feasibility alone never produces an output certificate.
    """
    if rounds < 1 or digits < 0:
        raise ValueError("rounds must be positive and digits nonnegative")
    lo, hi = np.array(problem.box, dtype=float).T
    start = (lo + hi) / 2 if seed is None else np.array(seed, dtype=float)
    if start.shape != (problem.model.num_params,) or not np.all(np.isfinite(start)):
        raise ValueError("seed must contain model.num_params finite parameters")
    xs = np.array(problem.sensors, dtype=float)
    y = np.array(problem.data, dtype=float)
    limit = 0.9 * float(problem.epsilon)
    k = problem.target
    others = [j for j in range(problem.model.num_params) if j != k]

    def values(t):
        return chebval(2 * xs - 1, solve_bvp(t, n=n, model=problem.model))

    def feasible(t):
        return (np.all(t >= lo) and np.all(t <= hi)
                and np.max(np.abs(values(t) - y)) <= limit)

    if not feasible(start):
        raise RuntimeError("seed is outside the prior or lacks the numerical safety margin")
    endpoints = []
    for sign in (-1, 1):
        t = start.copy()
        for _ in range(rounds):
            base = values(t)
            jac = np.empty((len(xs), problem.model.num_params))
            for j in range(problem.model.num_params):
                h = max(1e-3, abs(t[j]) * 1e-4)
                probe = t.copy()
                probe[j] += h
                jac[:, j] = (values(probe) - base) / h
            compensated = np.zeros(problem.model.num_params)
            compensated[k] = sign
            compensated[others] = np.linalg.lstsq(
                jac[:, others], -sign * jac[:, k], rcond=None
            )[0]
            direct = np.zeros(problem.model.num_params)
            direct[k] = sign
            for direction in (compensated, direct):
                caps = [(hi[j] - t[j]) / v if v > 0 else (lo[j] - t[j]) / v
                        for j, v in enumerate(direction) if v != 0]
                left, right = 0.0, max(0.0, min(caps))
                for _ in range(28):
                    step = (left + right) / 2
                    if feasible(t + step * direction):
                        left = step
                    else:
                        right = step
                t += left * direction
        theta = _rounded(t, digits)
        p = Poly(to_candidate(
            solve_bvp(theta, n=n, model=problem.model), boundary=problem.model.boundary
        ))
        report = check_witness(problem, theta, p, depth)
        if (report.verdict is not Verdict.VALID
                or Fraction(report.details['rho']) / 8 > problem.epsilon / 10):
            raise RuntimeError(f"search endpoint did not certify with residual budget: {report}")
        endpoints.append((theta, p))
    gap = abs(endpoints[0][0][k] - endpoints[1][0][k])
    if gap == 0:
        raise RuntimeError("search found no positive target gap")
    cert = make_pair_certificate(
        problem, *endpoints, gap,
        provenance={"generator": f"twindex.search {__version__}",
                    "notes": "Local heuristic; no global optimality claim."},
    )
    if check_certificate(cert, depth).verdict is not Verdict.VALID:
        raise RuntimeError("final pair certificate failed exact replay")
    return cert

sivia(problem, max_leaves=64, min_width='1/100', *, n=32, depth=6, candidate_factory=None, affine=None, mode='auto', priority='width', branching_policy='widest')

Build a bounded adaptive paving, then replay every claim exactly.

Bisect the widest coordinate (target first on ties), processing boxes with the widest target projection first. Numerical failures and budget-limited boxes remain unresolved, so even an unfinished search has an honest upper bound. candidate_factory(center) may supply an alternative exact Poly. mode selects auto (monotone when its hypotheses certify, else affine), monotone, affine, or plain. The legacy affine flag explicitly selects affine/plain and cannot be combined with a non-auto mode. Affine mode tries undecided plain boxes with centred finite-difference sensitivities of exact candidates. Probes stay in the prior; fixed coordinates use zero sensitivities. Failed probes or oversized expansions fall back to the plain result. No convergence or tightness is promised by this heuristic. priority="boundary" processes boxes nearest the prior target endpoints first, to focus on tightening the outer diameter.

Source code in twindex/search.py

def sivia(problem, max_leaves=64, min_width="1/100", *, n=32, depth=6,
          candidate_factory=None, affine=None, mode="auto", priority="width",
          branching_policy="widest"):
    """Build a bounded adaptive paving, then replay every claim exactly.

    Bisect the widest coordinate (target first on ties), processing boxes with
    the widest target projection first. Numerical failures and budget-limited
    boxes remain unresolved, so even an unfinished search has an honest upper
    bound. ``candidate_factory(center)`` may supply an alternative exact Poly.
    ``mode`` selects auto (monotone when its hypotheses certify, else affine),
    monotone, affine, or plain. The legacy ``affine`` flag explicitly selects
    affine/plain and cannot be combined with a non-auto mode.
    Affine mode tries undecided plain boxes with centred
    finite-difference sensitivities of exact candidates. Probes stay in the
    prior; fixed coordinates use zero sensitivities. Failed probes or oversized
    expansions fall back to the plain result. No convergence or tightness is
    promised by this heuristic. ``priority="boundary"`` processes boxes nearest
    the prior target endpoints first, to focus on tightening the outer diameter.
    """
    import heapq
    from itertools import count

    from .affine import sensor_range
    from .certificate import MAX_PAVING_DEPTH, MAX_PAVING_NODES, make_paving_certificate
    from .checker import check_box, monotone_hypotheses

    if branching_policy not in POLICIES:
        raise ValueError(f"branching_policy must be one of {POLICIES}")

    def split_forward(theta):
        coeffs = solve_bvp(theta, n=n, model=problem.model)
        return chebval(2*np.asarray(problem.sensors, dtype=float)-1, coeffs)

    if priority not in ("width", "boundary"):
        raise ValueError("priority must be width or boundary")

    if mode not in ("auto", "monotone", "affine", "plain"):
        raise ValueError("mode must be auto, monotone, affine or plain")
    if affine is not None:
        if not isinstance(affine, bool) or mode != "auto":
            raise ValueError("affine must be a bool and cannot be combined with mode")
        mode = "affine" if affine else "plain"
    if mode in ("auto", "monotone"):
        verdict, reason = monotone_hypotheses(problem, depth)
        if mode == "monotone" and verdict is not Verdict.VALID:
            raise ValueError(reason)
        mode = "monotone" if verdict is Verdict.VALID else "affine"
    width_limit = parse_rational(min_width)
    if (isinstance(max_leaves, bool) or not isinstance(max_leaves, int)
            or not 1 <= max_leaves <= (MAX_PAVING_NODES + 1) // 2):
        raise ValueError("max_leaves must be a positive integer within the schema node limit")
    if width_limit <= 0 or depth < 0:
        raise ValueError("min_width must be positive and depth nonnegative")
    if candidate_factory is None:
        def candidate_factory(center):
            return Poly(to_candidate(
                solve_bvp(center, n=n, model=problem.model), boundary=problem.model.boundary
            ))

    corner_cache = {}

    def corner_candidate(theta):
        if theta not in corner_cache:
            candidate = candidate_factory(theta)
            if not isinstance(candidate, Poly):
                raise ValueError("candidate_factory must return a Poly")
            corner_cache[theta] = candidate
        return corner_cache[theta]

    tree = {}
    queue = []
    serial = count()
    candidates = []
    indices = {}
    finished = []
    leaf_count = 1
    k = problem.target

    def enqueue(box, node, level):
        width = box[k][1] - box[k][0]
        distance = min(box[k][0] - problem.box[k][0], problem.box[k][1] - box[k][1])
        rank = (distance if priority == "boundary" else Fraction(0), -width)
        heapq.heappush(queue, (rank, next(serial), box, node, level))

    def classify(center, half, p, sens=(), corners=None):
        if sens:
            ranges = [sensor_range(p, sens, half, x) for x in problem.sensors]
            # If even zero residual cannot close either sufficient test, no
            # expansion can help. This shortcut only retains unresolved boxes.
            could_be_outside = any(lo - y > problem.epsilon or y - hi > problem.epsilon
                                   for (lo, hi), y in zip(ranges, problem.data, strict=True))
            could_be_inside = all(max(abs(lo - y), abs(hi - y)) <= problem.epsilon
                                  for (lo, hi), y in zip(ranges, problem.data, strict=True))
            if not could_be_outside and not could_be_inside:
                return "unresolved"
        report = check_box(problem, center, half, p, "inside", depth,
                           sensitivities=sens, corners=corners)
        if report.verdict is Verdict.VALID:
            return "inside"
        # INCONCLUSIVE for inside also means outside is undecided. Only
        # retry INVALID reports: these may prove disjointness or reject
        # malformed input. Never interpret INVALID alone as exclusion.
        if (report.verdict is Verdict.INVALID
                and check_box(problem, center, half, p, "outside", depth,
                              sensitivities=sens, corners=corners).verdict is Verdict.VALID):
            return "outside"
        return "unresolved"

    enqueue(problem.box, tree, 0)
    while queue:
        _, _, box, node, level = heapq.heappop(queue)
        center = tuple((lo + hi) / 2 for lo, hi in box)
        half = tuple((hi - lo) / 2 for lo, hi in box)
        status = "unresolved"
        sensitivities = ()
        corners = None

        failures = (ValueError, RuntimeError, ArithmeticError, np.linalg.LinAlgError)
        try:
            if mode == "monotone":
                corners = (corner_candidate(tuple(lo for lo, _ in box)),
                           corner_candidate(tuple(hi for _, hi in box)))
                p = None
                # Numerical screen can only retain a box as unresolved. Every
                # inside/outside label still requires a VALID exact report.
                xs = np.array(problem.sensors, dtype=float)
                y = np.array(problem.data, dtype=float)
                lower = np.polynomial.polynomial.polyval(
                    xs, [float(v) for v in corners[1].coeffs] or [0.0])
                upper = np.polynomial.polynomial.polyval(
                    xs, [float(v) for v in corners[0].coeffs] or [0.0])
                eps = float(problem.epsilon)
                margin = eps / 10
                likely_outside = np.any(lower > y + eps - margin) or np.any(
                    upper < y - eps + margin)
                likely_inside = np.all(lower >= y - eps - margin) and np.all(
                    upper <= y + eps + margin)
                if likely_outside or likely_inside:
                    wanted = "outside" if likely_outside else "inside"
                    report = check_box(problem, center, half, None, wanted, depth,
                                       corners=corners)
                    if report.verdict is Verdict.VALID:
                        status = wanted
                    elif report.verdict is Verdict.INVALID:
                        status = classify(center, half, None, corners=corners)
            else:
                p = candidate_factory(center)
                if not isinstance(p, Poly):
                    raise ValueError("candidate_factory must return a Poly")
                status = classify(center, half, p)
            if status == "unresolved" and mode == "affine":
                sens = []
                for j, h in enumerate(half):
                    if h == 0:
                        sens.append(Poly())
                        continue
                    step = min(Fraction(1, 10), h)
                    plus, minus = list(center), list(center)
                    plus[j] += step
                    minus[j] -= step
                    sp, sm = candidate_factory(tuple(plus)), candidate_factory(tuple(minus))
                    if not isinstance(sp, Poly) or not isinstance(sm, Poly):
                        raise ValueError("candidate_factory must return a Poly")
                    sens.append((sp - sm).scale(1 / (2 * step)))
                status = classify(center, half, p, sens)
                if status != "unresolved":
                    sensitivities = tuple(sens)
        except failures:
            # Failed numerics or an oversized expansion never exclude a box.
            pass
        if status != "unresolved":
            def register(poly):
                key = poly.coeffs
                if key not in indices:
                    indices[key] = len(candidates)
                    candidates.append(poly)
                return indices[key]

            if corners is not None:
                node.update(status=status, corners=[register(poly) for poly in corners])
            else:
                node.update(status=status, p=register(p))
            if sensitivities:
                node["s"] = [register(s) for s in sensitivities]
        else:
            widths = [hi - lo for lo, hi in box]
            axis = choose_axis(box, problem.box, k, branching_policy,
                               forward=split_forward, noise=[problem.epsilon]*len(problem.sensors),
                               target_ties=True)
            if (axis is not None and leaf_count < max_leaves and widths[axis] > width_limit
                    and level < MAX_PAVING_DEPTH):
                midpoint = center[axis]
                children = [{}, {}]
                node.update(split=axis, at=midpoint, children=children)
                for side, child in enumerate(children):
                    child_box = list(box)
                    lo, hi = box[axis]
                    child_box[axis] = (lo, midpoint) if side == 0 else (midpoint, hi)
                    enqueue(tuple(child_box), child, level + 1)
                leaf_count += 1
                continue
            node.update(status="unresolved")
        finished.append((box, status))

    def diameter(boxes):
        return max(b[k][1] for b in boxes) - min(b[k][0] for b in boxes)

    inside = [box for box, status in finished if status == "inside"]
    outer = [box for box, status in finished if status != "outside"]
    cert = make_paving_certificate(
        problem, candidates, tree, diameter(outer) if outer else Fraction(0),
        diameter(inside) if inside else None,
        provenance={"generator": f"twindex.search.sivia {__version__}",
                    "notes": f"Mode {mode}; priority {priority}; policy {branching_policy}; "
                             "unresolved leaves retained."},
    )
    report = check_certificate(cert, depth)
    if report.verdict is not Verdict.VALID:
        raise RuntimeError(f"final paving failed exact replay: {report.reasons}")
    return cert

twindex.solver

Untrusted floating-point Chebyshev collocation for the model BVP.

Coefficients use numpy's convention: sum c[k] T_k(2*x-1). Numerical convergence is only a diagnostic; only the exact checker can certify a solution.

residual(coeffs, theta, x, model=V1_MODEL)

Evaluate -p'' + q*f(p) - g at physical points x (untrusted).

Source code in twindex/solver.py

def residual(coeffs, theta, x, model=V1_MODEL):
    """Evaluate -p'' + q*f(p) - g at physical points x (untrusted)."""
    theta = _theta(theta, model)
    x = np.asarray(x, dtype=float)
    t = 2 * x - 1
    return (-cheb.chebval(t, cheb.chebder(coeffs, 2, scl=2))
            + _q(theta, x, model) * _f(cheb.chebval(t, coeffs), model)
            - _source(x, model))

solve_bvp(theta, n=48, tol=1e-13, max_iter=50, *, model=V1_MODEL)

Return n+1 Chebyshev coefficients, using n+1 Lobatto nodes.

Newton's equations include both boundary conditions. A backtracking line search requires a decreasing equation residual. Raise RuntimeError when the iteration budget is exhausted; low degree can still underresolve the continuous equation, so check the residual off the collocation grid too.

Source code in twindex/solver.py

def solve_bvp(theta, n=48, tol=1e-13, max_iter=50, *, model=V1_MODEL):
    """Return n+1 Chebyshev coefficients, using n+1 Lobatto nodes.

    Newton's equations include both boundary conditions. A backtracking line
    search requires a decreasing equation residual. Raise RuntimeError when
    the iteration budget is exhausted; low degree can still underresolve the
    continuous equation, so check the residual off the collocation grid too.
    """
    theta = _theta(theta, model)
    n, max_iter = operator.index(n), operator.index(max_iter)
    if n < 2 or max_iter < 1 or not np.isfinite(tol) or tol <= 0:
        raise ValueError("require n >= 2, max_iter >= 1 and finite tol > 0")
    t = np.cos(np.pi * np.arange(n + 1) / n)
    x = (t + 1) / 2
    v = cheb.chebvander(t, n)
    d2 = np.column_stack([
        cheb.chebval(t, cheb.chebder(np.eye(n + 1)[k], 2, scl=2))
        for k in range(n + 1)
    ])
    q = _q(theta, x, model)
    g = _source(x, model)
    c = np.zeros(n + 1)
    alpha, beta = map(float, model.boundary)
    # Linear boundary lift plus the Poisson solution for the midpoint source.
    midpoint_source = float(_source(np.array(0.5), model))
    c[0] = (alpha + beta) / 2 + midpoint_source / 16
    c[1] = (beta - alpha) / 2
    c[2] = -midpoint_source / 16

    def equations(c):
        u = v @ c
        r = -d2 @ c + q * _f(u, model) - g
        r[0], r[-1] = u[0] - beta, u[-1] - alpha
        return r

    for _ in range(max_iter):
        r = equations(c)
        norm = np.max(np.abs(r))
        if norm <= tol:
            return c
        u = v @ c
        jac = -d2 + (q * _f(u, model, derivative=True))[:, None] * v
        jac[0], jac[-1] = v[0], v[-1]
        try:
            step = np.linalg.solve(jac, -r)
        except np.linalg.LinAlgError as exc:
            raise RuntimeError("singular Newton system") from exc
        for backtrack in range(30):
            trial = c + 2.0**(-backtrack) * step
            if np.all(np.isfinite(trial)) and np.max(np.abs(equations(trial))) < norm:
                c = trial
                break
        else:
            raise RuntimeError(f"Newton line search stalled (residual {norm:.3g})")
    if np.max(np.abs(equations(c))) <= tol:
        return c
    raise RuntimeError(f"Newton did not converge within {max_iter} iterations")

twindex.candidate

Untrusted conversion to exact monomial candidates with exact prescribed BCs.

to_candidate(coeffs, bits=60, degree=None, *, boundary=(0, 0))

Return exact monomial coefficients, constant first.

Subtract the linear lift for the exact rational boundary=(alpha, beta). Divide the remaining floating Chebyshev polynomial by x(1-x), discard its boundary remainder, and round the quotient's Chebyshev coefficients to dyadics. Convert the quotient EXACTLY using T_(k+1)=2(2x-1)T_k-T_(k-1), then multiply by x(1-x). No floating monomial conversion is performed. degree optionally truncates the quotient before exact conversion. Add the lift back exactly, so endpoint values never depend on rounding.

Source code in twindex/candidate.py

def to_candidate(coeffs, bits=60, degree=None, *, boundary=(0, 0)):
    """Return exact monomial coefficients, constant first.

    Subtract the linear lift for the exact rational ``boundary=(alpha, beta)``.
    Divide the remaining floating Chebyshev polynomial by x(1-x), discard its boundary
    remainder, and round the quotient's Chebyshev coefficients to dyadics.
    Convert the quotient EXACTLY using T_(k+1)=2(2x-1)T_k-T_(k-1), then
    multiply by x(1-x). No floating monomial conversion is performed.
    ``degree`` optionally truncates the quotient before exact conversion.
    Add the lift back exactly, so endpoint values never depend on rounding.
    """
    bits = operator.index(bits)
    if bits < 0 or bits > 1024:
        raise ValueError("bits must lie between 0 and 1024")
    c = np.asarray(coeffs, dtype=float)
    if c.ndim != 1 or c.size == 0 or not np.all(np.isfinite(c)):
        raise ValueError("coeffs must be a nonempty finite one-dimensional array")
    alpha, beta = (parse_rational(v) for v in boundary)
    c = cheb.chebsub(c, [float((alpha + beta) / 2), float((beta - alpha) / 2)])
    quotient, _ = cheb.chebdiv(c, [1 / 8, 0, -1 / 8])
    if degree is not None:
        degree = operator.index(degree)
        if degree < 0:
            raise ValueError("degree must be nonnegative")
        quotient = quotient[:degree + 1]
    scale = 1 << bits
    # Fraction(float) is intentional here, outside the trusted core. Integer
    # rounding avoids overflow even when bits is larger than the float range.
    rounded = [Fraction(round(Fraction(float(a)) * scale), scale) for a in quotient]
    previous, current = [Fraction(1)], [Fraction(-1), Fraction(2)]
    s = [rounded[0]]
    for k in range(1, len(rounded)):
        s = _add(s, [rounded[k] * a for a in current])
        previous, current = current, _add(
            [2 * a for a in _shifted_t_times(current)], [-a for a in previous]
        )
    out = [Fraction(0)] * (len(s) + 2)
    for i, a in enumerate(s):
        out[i + 1] += a
        out[i + 2] -= a
    out[0] += alpha
    out[1] += beta - alpha
    while len(out) > 1 and out[-1] == 0:
        out.pop()
    return out

General numerical audit (untrusted)

These reports contain numerical evidence, not checker verdicts or certificates. See Audit any forward model for interpretation and examples.

twindex.general

Numerical evidence of ambiguity for user-supplied forward maps (UNTRUSTED).

These reports are not certificates and never carry checker verdicts. A small observed range, full-rank Jacobian, or unsuccessful search proves no uniqueness.

AuditReport dataclass

JSON-ready numerical evidence; endpoints have fresh forward evaluations.

parameters[k].observed_gap is a gap between sampled feasible points, never a global upper bound. Local LP ranges describe a linear approximation.

Source code in twindex/general.py

@dataclass
class AuditReport:
    """JSON-ready numerical evidence; endpoints have fresh forward evaluations.

    ``parameters[k].observed_gap`` is a gap between sampled feasible points,
    never a global upper bound. Local LP ranges describe a linear approximation.
    """

    problem: dict
    local: dict
    parameters: list[dict]
    search: dict
    evidence: str = "numerical_only"
    certified: bool = False
    limitations: tuple[str, ...] = (
        "Not a proof: forward error and floating-point error are not bounded.",
        "Search can miss disconnected feasible components and global extrema.",
        "Small observed gaps and full local rank do not prove identifiability.",
        "No feasible sample does not prove that the compatible set is empty.",
    )

    def typed(self) -> NumericalEvidence:
        return NumericalEvidence(self.to_dict())

    def to_dict(self):
        """Return a JSON-serializable report, distinct from certificate JSON."""
        return {"format": "twindex-numerical-audit", "version": 1, **asdict(self)}

to_dict()

Return a JSON-serializable report, distinct from certificate JSON.

Source code in twindex/general.py

def to_dict(self):
    """Return a JSON-serializable report, distinct from certificate JSON."""
    return {"format": "twindex-numerical-audit", "version": 1, **asdict(self)}

CertifiedAuditResult dataclass

Separate numerical snapshot and exact replay attempts, including failures.

Source code in twindex/general.py

@dataclass
class CertifiedAuditResult:
    """Separate numerical snapshot and exact replay attempts, including failures."""

    numerical: dict
    certified: list[dict]

    def typed(self, *, depth: int = 6) -> tuple[CertifiedResult | Inconclusive, ...]:
        """One independently replayed outcome per audited coordinate.

        The numerical snapshot stays separate at .numerical; no endpoint or
        failed generation is silently promoted to a certificate.
        """
        outcomes = []
        for section in self.certified:
            attempts = section.get('attempts', [])
            candidates = [a['certificate'] for a in attempts if 'certificate' in a]
            if candidates:
                outcomes.append(typed_certificate(candidates[-1], depth=depth))
            else:
                reasons = tuple(a['generation_error'] for a in attempts if 'generation_error' in a)
                outcomes.append(Inconclusive('INCONCLUSIVE_INSUFFICIENT_BOUND', reasons or
                                             ('No candidate certificate generated.',)))
        return tuple(outcomes)

    def to_dict(self):
        return {"format": "twindex-certified-audit", "version": 1,
                **asdict(self)}

typed(*, depth=6)

One independently replayed outcome per audited coordinate.

The numerical snapshot stays separate at .numerical; no endpoint or failed generation is silently promoted to a certificate.

Source code in twindex/general.py

def typed(self, *, depth: int = 6) -> tuple[CertifiedResult | Inconclusive, ...]:
    """One independently replayed outcome per audited coordinate.

    The numerical snapshot stays separate at .numerical; no endpoint or
    failed generation is silently promoted to a certificate.
    """
    outcomes = []
    for section in self.certified:
        attempts = section.get('attempts', [])
        candidates = [a['certificate'] for a in attempts if 'certificate' in a]
        if candidates:
            outcomes.append(typed_certificate(candidates[-1], depth=depth))
        else:
            reasons = tuple(a['generation_error'] for a in attempts if 'generation_error' in a)
            outcomes.append(Inconclusive('INCONCLUSIVE_INSUFFICIENT_BOUND', reasons or
                                         ('No candidate certificate generated.',)))
    return tuple(outcomes)

audit_model(forward, data, prior, epsilon, *, estimate=None, jacobian=None, starts=4, seed=0, max_iter=200, diff_step=1e-05, rank_rtol=None)

Audit any deterministic forward(theta) -> measurements callable.

Requires twindex[general]. data is a nonempty vector, prior an (N, 2) finite box (fixed coordinates allowed), and epsilon a positive scalar or per-measurement vector. Only bounded componentwise noise is used. A supplied jacobian(theta) has shape (M, N), in physical parameter units. Default relative rank tolerance is max(1e-7, diff_step) for finite differences, or max(M, Nfree)*machine_epsilon for supplied derivatives. rank_rtol overrides it; all rank diagnostics remain numerical evidence.

Parameters are scaled to their prior widths and residuals to epsilon. The reference point is the best bounded least-squares fit from starts starts (including estimate or the midpoint). Local SVD and LP diagnostics are followed by multistart nonlinear constrained endpoint searches. All accepted endpoints are re-evaluated without relaxing the noise or prior constraints. Solver success alone never establishes feasibility. Forward failures or malformed outputs raise; stochastic and stateful maps are unsupported.

Source code in twindex/general.py

def audit_model(forward, data, prior, epsilon, *, estimate=None, jacobian=None,
                starts=4, seed=0, max_iter=200, diff_step=1e-5, rank_rtol=None):
    """Audit any deterministic ``forward(theta) -> measurements`` callable.

    Requires ``twindex[general]``. ``data`` is a nonempty vector, ``prior`` an
    (N, 2) finite box (fixed coordinates allowed), and ``epsilon`` a positive
    scalar or per-measurement vector. Only bounded componentwise noise is used.
    A supplied ``jacobian(theta)`` has shape (M, N), in physical parameter units.
    Default relative rank tolerance is max(1e-7, diff_step) for finite differences,
    or max(M, Nfree)*machine_epsilon for supplied derivatives. ``rank_rtol``
    overrides it; all rank diagnostics remain numerical evidence.

    Parameters are scaled to their prior widths and residuals to epsilon. The
    reference point is the best bounded least-squares fit from ``starts`` starts
    (including ``estimate`` or the midpoint). Local SVD and LP diagnostics are
    followed by multistart nonlinear constrained endpoint searches. All accepted
    endpoints are re-evaluated without relaxing the noise or prior constraints.
    Solver success alone never establishes feasibility. Forward failures or
    malformed outputs raise; stochastic and stateful maps are unsupported.
    """
    from scipy.optimize import least_squares, linprog, minimize

    if not callable(forward) or (jacobian is not None and not callable(jacobian)):
        raise ValueError("forward and optional jacobian must be callable")
    y = _array(data, "data")
    if y.ndim != 1 or not y.size:
        raise ValueError("data must be a nonempty vector")
    box = _array(prior, "prior")
    if box.ndim != 2 or box.shape[1] != 2 or not len(box):
        raise ValueError("prior must have shape (N, 2), N >= 1")
    lo, hi = box.T
    with np.errstate(over="ignore"):
        width = hi - lo
    if np.any(width < 0) or not np.all(np.isfinite(width)):
        raise ValueError("prior widths must be finite and nonnegative")
    eps = _array(epsilon, "epsilon")
    if eps.ndim == 0:
        eps = np.full(y.shape, float(eps))
    if eps.shape != y.shape or np.any(eps <= 0):
        raise ValueError("epsilon must be positive and scalar or match data")
    for value, name in ((starts, "starts"), (max_iter, "max_iter")):
        if isinstance(value, bool) or not isinstance(value, int) or value < 1:
            raise ValueError(f"{name} must be a positive integer")
    if isinstance(seed, bool) or not isinstance(seed, int) or seed < 0:
        raise ValueError("seed must be a nonnegative integer")
    if not np.isfinite(diff_step) or not 0 < diff_step < 0.5:
        raise ValueError("diff_step must be finite and in (0, 0.5)")
    if rank_rtol is not None and (not np.isfinite(rank_rtol) or not 0 < rank_rtol < 1):
        raise ValueError("rank_rtol must be finite and in (0, 1)")
    active = np.flatnonzero(width > 0)
    n = len(lo)
    nfree = len(active)
    theta0 = lo + width / 2 if estimate is None else _array(estimate, "estimate", (n,))
    if np.any(theta0 < lo) or np.any(theta0 > hi):
        raise ValueError("estimate must lie in prior")
    calls = 0
    runs = []
    samples = []

    def unpack(z):
        theta = lo.copy()
        theta[active] += width[active] * z
        return theta

    def evaluate(theta):
        nonlocal calls
        calls += 1
        # User solvers may reuse an output buffer on the next evaluation.
        pred = _array(forward(theta.copy()), "forward output", y.shape).copy()
        with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
            residual = (pred - y) / eps
        if not np.all(np.isfinite(residual)):
            raise ValueError("scaled forward residual must be finite")
        return pred, residual

    def residual(z):
        return evaluate(unpack(z))[1]

    def accept(z):
        # No clipping or reuse of optimizer residuals: audit exactly this point.
        z = _array(z, "optimizer point", (nfree,))
        theta = unpack(z)
        if np.any(z < 0) or np.any(z > 1) or np.any(theta < lo) or np.any(theta > hi):
            return
        pred, r = evaluate(theta)
        misfit = float(np.max(np.abs(r)))
        if misfit <= 1:
            samples.append({"theta": theta.tolist(), "prediction": pred.tolist(),
                            "max_scaled_misfit": misfit})

    rng = np.random.default_rng(seed)
    z0 = (theta0[active] - lo[active]) / width[active]
    initial = [z0] + [rng.random(nfree) for _ in range(starts - 1)]
    fitted = []
    for z in initial:
        accept(z)
        if nfree:
            fit = least_squares(residual, z, bounds=(0, 1), max_nfev=max_iter,
                                ftol=1e-10, xtol=1e-10, gtol=1e-10)
            z = fit.x
            runs.append({"stage": "fit", "success": bool(fit.success),
                         "message": str(fit.message)})
            accept(z)
        fitted.append(z)
    reference = min(fitted, key=lambda z: np.linalg.norm(residual(z)))
    theta_ref = unpack(reference)
    pred_ref, r_ref = evaluate(theta_ref)

    if jacobian is not None:
        j = _array(jacobian(theta_ref.copy()), "jacobian", (len(y), n))
        j = j[:, active] * width[active] / eps[:, None]
    else:
        j = np.zeros((len(y), nfree))
        for k in range(nfree):
            left, right = reference.copy(), reference.copy()
            left[k] = max(0, left[k] - diff_step)
            right[k] = min(1, right[k] + diff_step)
            if right[k] == left[k]:
                raise ValueError("finite difference step is unrepresentable")
            j[:, k] = (residual(right) - residual(left)) / (right[k] - left[k])
    if not np.all(np.isfinite(j)):
        raise ValueError("scaled jacobian must be finite")
    # Preserve all right-null directions when N > M without allocating an MxM
    # left factor for observation-heavy models (e.g. sampled solution fields).
    _, singular, vh = np.linalg.svd(j, full_matrices=nfree > len(y))
    if not np.all(np.isfinite(singular)):
        raise ValueError("singular values overflowed; rescale your model")
    rtol = (float(rank_rtol) if rank_rtol is not None else
            max(j.shape) * np.finfo(float).eps if jacobian is not None else
            max(1e-7, diff_step))
    ratios = singular / singular[0] if singular.size and singular[0] else np.zeros_like(singular)
    threshold = float(rtol * singular[0]) if singular.size else 0
    rank = int(np.count_nonzero(singular > threshold))
    null = np.zeros((nfree - rank, n))
    null[:, active] = vh[rank:]
    local = {"reference": theta_ref.tolist(), "prediction": pred_ref.tolist(),
             "max_scaled_misfit": float(np.max(np.abs(r_ref))),
             "active_parameters": active.tolist(), "scaled_jacobian": j.tolist(),
             "singular_values": singular.tolist(), "rank": rank,
             "rank_threshold": threshold, "rank_rtol": rtol,
             "derivative_source": "supplied" if jacobian is not None else "finite_difference",
             "singular_value_ratios": ratios.tolist(),
             "null_direction_meaning": "numerically rank-deficient directions (ratio <= rank_rtol)",
             "null_directions_scaled": null.tolist(),
             "linearised_ranges": []}

    # LP in d = z - reference, with dimensionless residual inequalities.
    a = np.vstack((j, -j))
    b = np.concatenate((1 - r_ref, 1 + r_ref))
    bounds = list(zip(-reference, 1 - reference, strict=True))
    lp_seeds = []
    for k in range(n):
        if width[k] == 0:
            local["linearised_ranges"].append({"parameter": k, "status": "fixed_by_prior",
                                                "range": [float(lo[k]), float(hi[k])]})
            continue
        index = int(np.flatnonzero(active == k)[0])
        c = np.zeros(nfree)
        c[index] = 1
        endpoints = []
        violations = []
        statuses = []
        for sign in (1, -1):
            result = linprog(sign * c, A_ub=a, b_ub=b, bounds=bounds, method="highs",
                             options={"primal_feasibility_tolerance": 1e-10,
                                      "dual_feasibility_tolerance": 1e-10})
            statuses.append(str(result.message))
            if result.success:
                d = result.x
                violation = max(0., float(np.max(a @ d - b)),
                                float(np.max(-reference - d)),
                                float(np.max(d - (1 - reference))))
                if violation <= 1e-8:
                    endpoints.append(float(unpack(reference + d)[k]))
                    violations.append(violation)
                    lp_seeds.append(np.clip(reference + d, 0, 1))
        local["linearised_ranges"].append({
            "parameter": k, "status": "solved" if len(endpoints) == 2 else "unavailable",
            "range": endpoints if len(endpoints) == 2 else None, "messages": statuses,
            "max_scaled_violation": max(violations) if violations else None,
        })

    # Both LP seeds and multistart fits are tried even when no feasible seed was
    # found. Local linear feasibility is never assumed to imply actual feasibility.
    seeds = fitted + lp_seeds
    for z in lp_seeds:
        accept(z)

    def constraints(z):
        r = residual(z)
        return np.concatenate((1 - r, 1 + r))

    for index, k in enumerate(active):
        for sign in (1, -1):
            objective_jac = np.zeros(nfree)
            objective_jac[index] = sign
            for z in seeds:
                result = minimize(lambda v, s=sign, i=index: s * v[i], z,
                                  jac=lambda v, c=objective_jac: c,
                                  method="SLSQP", bounds=[(0, 1)] * nfree,
                                  constraints={"type": "ineq", "fun": constraints},
                                  options={"maxiter": max_iter, "ftol": 1e-10})
                runs.append({"stage": "extreme", "parameter": int(k), "sign": sign,
                             "success": bool(result.success), "message": str(result.message)})
                accept(result.x)
    parameters = []
    for k in range(n):
        if samples:
            low = min(samples, key=lambda v: v["theta"][k])
            high = max(samples, key=lambda v: v["theta"][k])
            gap = high["theta"][k] - low["theta"][k]
            item = {"parameter": k, "status": "separated_samples_found" if gap > 0
                    else "no_separated_samples_found", "low": low, "high": high,
                    "observed_gap": gap,
                    "prior_fraction": gap / width[k] if width[k] else None}
        else:
            item = {"parameter": k, "status": "no_feasible_sample_found", "low": None,
                    "high": None, "observed_gap": None, "prior_fraction": None}
        parameters.append(item)
    return AuditReport(
        problem={"data": y.tolist(), "epsilon": eps.tolist(), "prior": box.tolist(),
                 "noise_model": "componentwise_bounded"},
        local=local, parameters=parameters,
        search={"starts": starts, "seed": seed, "max_iter": max_iter, "diff_step": diff_step,
                "forward_calls": calls, "feasible_samples": len(samples), "runs": runs},
    )

certify_audit(report, problem, theta_map=None, *, target_map=None, n=32, depth=6, shrink_factors=('1', '9/10', '1/2'), digits=12)

Replay numerical endpoints against an explicitly supplied exact 1-D Problem.

The caller asserts correspondence of the callable and exact problem. This function cannot establish that assertion. By default coordinates correspond in order. A custom theta_map requires target_map[k], the exact target index to certify for audit coordinate k. A map need not preserve separation: each claim uses the actual mapped rational gap. No upper bound is attempted.

Source code in twindex/general.py

def certify_audit(report, problem, theta_map=None, *, target_map=None, n=32,
                  depth=6, shrink_factors=("1", "9/10", "1/2"), digits=12):
    """Replay numerical endpoints against an explicitly supplied exact 1-D Problem.

    The caller asserts correspondence of the callable and exact problem. This
    function cannot establish that assertion. By default coordinates correspond
    in order. A custom theta_map requires target_map[k], the exact target index
    to certify for audit coordinate k. A map need not preserve separation: each
    claim uses the actual mapped rational gap. No upper bound is attempted.
    """
    from dataclasses import replace
    from fractions import Fraction

    from .candidate import to_candidate
    from .certificate import make_pair_certificate
    from .checker import Verdict, check_certificate
    from .model import Problem
    from .poly import Poly
    from .rational import parse_rational
    from .solver import solve_bvp

    if not isinstance(report, AuditReport) or not isinstance(problem, Problem):
        raise ValueError("requires an AuditReport and an explicit exact 1-D Problem")
    if isinstance(digits, bool) or not isinstance(digits, int) or not 1 <= digits <= 30:
        raise ValueError("digits must be an integer in 1..30")
    count = len(report.parameters)
    if theta_map is None:
        if count != problem.model.num_params:
            raise ValueError("coordinate counts differ; provide theta_map and target_map")
        def theta_map(theta):
            return theta
        if target_map is None:
            target_map = list(range(count))
    if not callable(theta_map) or target_map is None or len(target_map) != count:
        raise ValueError("custom theta_map requires one exact target_map index per parameter")
    problems = [replace(problem, target=k) for k in target_map]
    factors = tuple(parse_rational(f) for f in shrink_factors)
    if not factors or any(not 0 < f <= 1 for f in factors):
        raise ValueError("shrink_factors must be exact rationals in (0, 1]")

    def rationalise(theta):
        values = _array(theta_map(np.array(theta, dtype=float).copy()), "theta_map output",
                        (problem.model.num_params,))
        return tuple(Fraction(format(float(v), f".{digits}f")) for v in values)

    reference = rationalise(report.local["reference"])
    sections = []
    for item, exact_problem in zip(report.parameters, problems, strict=True):
        section = {"audit_parameter": item["parameter"], "exact_target": exact_problem.target,
                   "attempts": [], "status": "no_endpoints"}
        sections.append(section)
        if item["low"] is None or item["high"] is None:
            continue
        endpoints = [rationalise(item[key]["theta"]) for key in ("low", "high")]
        section["status"] = "not_certified"
        for factor in factors:
            thetas = [tuple(r + factor * (v - r)
                            for r, v in zip(reference, endpoint, strict=True))
                      for endpoint in endpoints]
            attempt = {"shrink_factor": str(factor),
                       "theta": [[str(v) for v in theta] for theta in thetas]}
            section["attempts"].append(attempt)
            try:
                witnesses = []
                for theta in thetas:
                    p = Poly(to_candidate(solve_bvp(theta, n=n, model=problem.model),
                                          boundary=problem.model.boundary))
                    witnesses.append((theta, p))
                k = exact_problem.target
                gap = abs(witnesses[1][0][k] - witnesses[0][0][k])
                cert = make_pair_certificate(exact_problem, *witnesses, gap)
                replay = check_certificate(cert, depth=depth)
                attempt.update(certificate=cert, replay=replay.to_dict())
                if replay.verdict is Verdict.VALID:
                    section["status"] = "certified_pair"
                    break
            except (ValueError, RuntimeError, ArithmeticError, np.linalg.LinAlgError) as error:
                attempt["generation_error"] = str(error)
    return CertifiedAuditResult(report.to_dict(), sections)