"""The one claim that survives the objection a reader raises first. Everything else in this project is a tie, a loss, or a negative result. This file pins the exception and, just as important, pins its scope: it is a claim about REPAIRING regula falsi, not about being the best root finder. """ from __future__ import annotations import math import os import sys import pytest sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from simroot.baselines import (anderson_bjorck, brent, chandrupatla, illinois, pegasus, regula_falsi) from stage18_bulletproof import falsi_driver, load_roots XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6 PROBS = load_roots() REPAIRS = ("none", "illinois", "pegasus", "anderson", "credit") def score(method): total = solved = worst = 0 for f, a, b, root in PROBS: try: x, n, conv = method(f, a, b, XTOL, MAXFEV) except Exception: x, n, conv = float("nan"), MAXFEV, False ok = bool(conv and math.isfinite(x) and abs(x - root) <= ACC) total += n if ok else MAXFEV solved += ok if ok: worst = max(worst, n) return total, solved, worst def by_repair(repair): return score(lambda f, a, b, xt, mx: falsi_driver(f, a, b, xt, mx, repair)) # -------------------------------------------------------------------------- # the driver cannot be handicapping the competitors # -------------------------------------------------------------------------- @pytest.mark.parametrize("repair,independent", [ ("none", regula_falsi), ("illinois", illinois), ("pegasus", pegasus), ("anderson", anderson_bjorck)]) def test_two_independent_implementations_agree_exactly(repair, independent): """The classical repairs exist twice in this repository, written at different times against different structures: `_false_position_family` in simroot.baselines, and the pluggable driver in stage18. They agree to the evaluation, which is what removes 'your driver favours your repair' from the table.""" a = score(independent) b = by_repair(repair) assert a[:2] == b[:2], (repair, a, b) # -------------------------------------------------------------------------- # the claim # -------------------------------------------------------------------------- def test_credit_beats_every_classical_repair_on_both_axes_at_once(): """Same interpolation, same driver, same budget; only the repair differs. Against Illinois (1971), Pegasus (1972) and Anderson-Bjorck (1973).""" c_total, c_solved, c_worst = by_repair("credit") for repair in ("illinois", "pegasus", "anderson"): t, s, w = by_repair(repair) assert c_solved > s, (repair, c_solved, s) assert c_worst < w, (repair, c_worst, w) assert c_total < min(by_repair(r)[0] for r in ("illinois", "pegasus", "anderson")) def test_no_classical_repair_gets_near_the_solved_count(): """193 of 199 is what bisection, Chandrupatla and TOMS 748 achieve; the best classical repair of regula falsi reaches 176.""" _t, c_solved, _w = by_repair("credit") assert c_solved == 193, c_solved for repair in ("none", "illinois", "pegasus", "anderson"): assert by_repair(repair)[1] <= 176, repair def test_the_advantage_is_not_an_artefact_of_aggregating_two_sets(): """A reader will ask whether the win comes from one set. On the wide-bracket set the credit repair halves the cost and doubles the solved count; on the canonical set it matches Illinois's solved count at lower cost. Anderson-Bjorck is cheaper there, and fails eight problems the credit repair solves - that trade is stated, not hidden.""" from scipy.optimize._tstutils import (_APS_TESTS_DICTS, _CHANDRUPATLA_TESTS_DICTS) def per_set(repair, dicts, chand): total = solved = 0 for d in dicts: a, b = float(d["bracket"][0]), float(d["bracket"][1]) args = () if chand else d.get("args", ()) f = (lambda ff, aa: lambda t: ff(t, *aa))(d["f"], args) try: x, n, conv = falsi_driver(f, a, b, XTOL, MAXFEV, repair) except Exception: x, n, conv = float("nan"), MAXFEV, False ok = bool(conv and math.isfinite(x) and abs(x - float(d["root"])) <= ACC) total += n if ok else MAXFEV solved += ok return total, solved for repair in ("illinois", "pegasus", "anderson"): c_t, c_s = per_set("credit", _CHANDRUPATLA_TESTS_DICTS, True) t, s = per_set(repair, _CHANDRUPATLA_TESTS_DICTS, True) assert c_t < 0.7 * t, (repair, c_t, t) # halves the cost assert c_s >= 1.7 * s, (repair, c_s, s) # doubles the solved def test_the_scope_of_the_claim_is_repairs_not_the_field(): """The credit repair of regula falsi is still beaten by Brent and by Chandrupatla. The claim must never be stated as 'the best root finder'.""" c_total, _s, _w = by_repair("credit") for better in (brent, chandrupatla): t, _s2, _w2 = score(lambda f, a, b, xt, mx, _m=better: _m(f, a, b, xt, mx)) assert t < c_total, (better.__name__, t, c_total) # -------------------------------------------------------------------------- # the claim that does NOT survive, pinned so it is not overstated # -------------------------------------------------------------------------- def test_the_anc_claim_is_not_uniform_and_must_be_stated_as_such(): """Against a Kwong-Johnston oracle tuned per scenario the measured step size wins on the mean (19.14 dB against 17.51) but on only four of ten cells, and the mean is carried by two of them. It is a real result and it is not bullet-proof; this test exists so the difference is not blurred.""" import numpy as np from anc.bench import SCENARIOS, make_primary_path, make_secondary_path, \ run_nlms from anc.rules import Measured from stage18_bulletproof import KJTuned, KJ_GRID wins = cells = 0 for _sname, gen in SCENARIOS: for snr in (30.0, 15.0): def sc(rule): vals = [] for s in range(3): rng = np.random.default_rng(1000 + s) noise = gen(8000, rng) nr = run_nlms(noise, make_primary_path(rng), make_secondary_path(rng), rule, snr_db=snr, seed=s) if math.isfinite(nr): vals.append(nr) return sum(vals) / len(vals) ours = sc(Measured()) best_kj = max(sc(KJTuned(a, g)) for a, g in KJ_GRID) cells += 1 wins += ours > best_kj assert 0 < wins < cells, (wins, cells)