vvsandClaude Opus 5
68df1d25db
Impact, every sweep regenerated: the paper's key configurations are unchanged to the evaluation (proposed 2806/3220, credit 3568/3020, Chandrupatla, TOMS 748, Brent, bisection); 46 stated numbers move (ITP CHA 2857 -> 2896, ITP oracle 2611 -> 2647, record CHA figures by 14-18, the grid means by ~1%) and one finding weakens: the measured-order exponent of section 7.1 was cheaper in 29 of 60 far cells and 7 of 7 safeguards; with the fix 12 of 61 and 6 of 7, and on the wide-bracket set it now loses in every safeguard and drops 5 problems. PREREG G1-G3 outcome on the fixed code: G2 (repeating step pattern, period from the data) PASSES Y1-Y4 and is ACCEPTED by the stop rule - proposed rule APS 2806 -> 2782 (below Chandrupatla's 2787), CHA 3220 -> 3169, both 6026 -> 5951, per problem vs Chandrupatla on APS 95/20/38 p = 8e-7. G1 fails Y1/Y2 (its Y3 miss was the sign bug), G3 fails Y2, A fails Y2. Regression 887 pass; the two red tests are the paper-number check (awaiting the paper update) and test_blend (fixed since, copy in stage16). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
231 lines
9.1 KiB
Python
231 lines
9.1 KiB
Python
"""Stage 18 - attack the two claims that look like wins.
|
|
|
|
A claim is worth publishing only if it survives the objection a reader will
|
|
raise first. For each of the two candidate wins the objection is obvious, and
|
|
this file runs it rather than answering it in prose.
|
|
|
|
CLAIM A. The credit safeguard rescues unsafe interpolation rules: regula
|
|
falsi goes from 133 solved and 18799 evaluations to 193 and 8056.
|
|
OBJECTION: nobody uses raw regula falsi. It has been fixed since 1971 -
|
|
Illinois, Pegasus, Anderson-Bjorck - and those fixes are the real
|
|
comparison. Beating the unfixed version proves nothing.
|
|
TEST: the same rule, the same driver, the same everything, with the
|
|
stale-endpoint correction swapped for each classical fix in turn and for
|
|
the credit projection. Only the repair differs.
|
|
|
|
CLAIM B. In active noise control the measured step size beats Kwong-Johnston
|
|
by 7 dB.
|
|
OBJECTION: Kwong-Johnston's gains are scale-dependent, and yours are the
|
|
published defaults run on a signal they were not scaled for. Give it its
|
|
best constants and the gap will close.
|
|
TEST: an oracle over its (alpha, gamma) grid, chosen separately for every
|
|
scenario - not an achievable controller, exactly as the ITP oracle in the
|
|
root-finding work was not.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
from simroot.baselines import _opposite, _same_sign # 2026-09-21 fix
|
|
|
|
import math
|
|
|
|
import numpy as np
|
|
from scipy.optimize._tstutils import (_APS_TESTS_DICTS,
|
|
_CHANDRUPATLA_TESTS_DICTS)
|
|
|
|
from anc.bench import SCENARIOS, make_primary_path, make_secondary_path, run_nlms
|
|
from anc.rules import Fixed, KwongJohnston, Measured, Rule
|
|
from simroot.baselines import Counter, _Budget
|
|
from simroot.pure import _Safeguard
|
|
|
|
XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
|
|
|
|
|
|
# ==========================================================================
|
|
# CLAIM A: the credit projection against the classical repairs of regula falsi
|
|
# ==========================================================================
|
|
def falsi_driver(f, a, b, xtol, maxfev, repair):
|
|
"""Regula falsi with a pluggable repair for the stale endpoint.
|
|
|
|
Everything is shared: the interpolation, the bracket update, the budget,
|
|
the duplicate handling. `repair` is the only difference, so a difference
|
|
in the result cannot come from anywhere else.
|
|
|
|
"none" the unfixed rule
|
|
"illinois" halve the stale ordinate (Hafner/Dowell 1971)
|
|
"pegasus" f * fx / (fx + f_new) (Dowell/Jarratt 1972)
|
|
"anderson" f * (1 - f_new/f_moved) (Anderson/Bjorck 1973)
|
|
"credit" leave the ordinates alone; project the step onto the
|
|
interval the earned bisection budget allows
|
|
"""
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b)
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
if _same_sign(fa, fb):
|
|
raise ValueError("no sign change")
|
|
sg = _Safeguard("credit", abs(b - a), xtol)
|
|
while abs(b - a) > xtol:
|
|
lo, hi = (a, b) if a < b else (b, a)
|
|
denom = fb - fa
|
|
x_new = 0.5 * (a + b) if denom == 0.0 else (a * fb - b * fa) / denom
|
|
if repair == "credit":
|
|
x_new = sg.project_credit(x_new, lo, hi)
|
|
if not math.isfinite(x_new) or not (lo < x_new < hi):
|
|
x_new = 0.5 * (a + b)
|
|
if x_new in (a, b):
|
|
x_new = 0.5 * (a + b)
|
|
if x_new in (a, b):
|
|
return x_new, g.n, abs(b - a) <= xtol
|
|
w_before = abs(b - a)
|
|
fx = g(x_new)
|
|
if fx == 0.0:
|
|
return x_new, g.n, True
|
|
if _opposite(fa, fx):
|
|
if repair == "illinois":
|
|
fa *= 0.5
|
|
elif repair == "pegasus":
|
|
fa = fa * fb / (fb + fx)
|
|
elif repair == "anderson":
|
|
m = 1.0 - fx / fb
|
|
fa = fa * (m if m > 0.0 else 0.5)
|
|
b, fb = x_new, fx
|
|
else:
|
|
if repair == "illinois":
|
|
fb *= 0.5
|
|
elif repair == "pegasus":
|
|
fb = fb * fa / (fa + fx)
|
|
elif repair == "anderson":
|
|
m = 1.0 - fx / fa
|
|
fb = fb * (m if m > 0.0 else 0.5)
|
|
a, fa = x_new, fx
|
|
if repair == "credit":
|
|
sg.note_shrink(w_before, abs(b - a))
|
|
x = x_new
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
def load_roots():
|
|
out = []
|
|
for d in _APS_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
args = d.get("args", ())
|
|
out.append(((lambda ff, aa: lambda t: ff(t, *aa))(d["f"], args),
|
|
float(a), float(b), float(d["root"])))
|
|
for d in _CHANDRUPATLA_TESTS_DICTS:
|
|
a, b = d["bracket"]
|
|
out.append(((lambda ff: lambda t: ff(t))(d["f"]),
|
|
float(a), float(b), float(d["root"])))
|
|
return out
|
|
|
|
|
|
def claim_a():
|
|
probs = load_roots()
|
|
print("=" * 84)
|
|
print("CLAIM A - the credit projection against the classical repairs")
|
|
print("=" * 84)
|
|
print("same interpolation, same driver, same budget; only the repair "
|
|
"differs\n")
|
|
print(f"{'repair':16s}{'year':>7s}{'evals':>10s}{'solved':>10s}"
|
|
f"{'worst':>8s}{'over bisect+3':>15s}")
|
|
print("-" * 84)
|
|
rows = {}
|
|
for repair, year in (("none", "300 BC"), ("illinois", "1971"),
|
|
("pegasus", "1972"), ("anderson", "1973"),
|
|
("credit", "here")):
|
|
total = solved = worst = over = 0
|
|
for f, a, b, root in probs:
|
|
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 - root) <= ACC)
|
|
total += n if ok else MAXFEV
|
|
solved += ok
|
|
if ok:
|
|
worst = max(worst, n)
|
|
if n > math.ceil(math.log2(abs(b - a) / XTOL)) + 4:
|
|
over += 1
|
|
rows[repair] = (total, solved, worst, over)
|
|
print(f"{repair:16s}{year:>7s}{total:>10d}{solved:>7d}/199"
|
|
f"{worst:>8d}{over:>15d}", flush=True)
|
|
best_classical = min(rows[r][0] for r in ("illinois", "pegasus",
|
|
"anderson"))
|
|
c = rows["credit"][0]
|
|
print(f"\n best classical repair {best_classical}, credit {c} "
|
|
f"=> {'SURVIVES' if c < best_classical else 'FAILS'} the objection")
|
|
return rows
|
|
|
|
|
|
# ==========================================================================
|
|
# CLAIM B: the measured step size against a Kwong-Johnston oracle
|
|
# ==========================================================================
|
|
class KJTuned(KwongJohnston):
|
|
def __init__(self, a, g):
|
|
self.A, self.G = a, g
|
|
self.name = f"KJ a={a:g} g={g:g}"
|
|
self.constants = 4
|
|
|
|
|
|
N, SEEDS, SNRS = 8000, 6, (30.0, 15.0)
|
|
KJ_GRID = [(a, g) for a in (0.9, 0.97, 0.99)
|
|
for g in (2e-4, 2e-3, 2e-2, 2e-1)]
|
|
|
|
|
|
def anc_score(rule, scen_gen, snr, seeds=SEEDS):
|
|
vals = []
|
|
for s in range(seeds):
|
|
rng = np.random.default_rng(1000 + s)
|
|
noise = scen_gen(N, 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)) if vals else float("nan")
|
|
|
|
|
|
def claim_b():
|
|
print("\n" + "=" * 84)
|
|
print("CLAIM B - the measured step size against a Kwong-Johnston ORACLE")
|
|
print("=" * 84)
|
|
print(f"the oracle is given the best of {len(KJ_GRID)} (alpha, gamma) "
|
|
f"pairs for every scenario\nand noise level separately - not an "
|
|
f"achievable controller\n")
|
|
print(f"{'scenario':22s}{'SNR':>6s}{'KJ default':>12s}{'KJ oracle':>12s}"
|
|
f"{'measured':>11s}{'fixed 0.5':>11s}")
|
|
print("-" * 84)
|
|
tot = {"kj": 0.0, "orc": 0.0, "meas": 0.0, "fix": 0.0}
|
|
cells = 0
|
|
for sname, gen in SCENARIOS:
|
|
for snr in SNRS:
|
|
kj = anc_score(KwongJohnston(), gen, snr)
|
|
orc = max(anc_score(KJTuned(a, g), gen, snr) for a, g in KJ_GRID)
|
|
meas = anc_score(Measured(), gen, snr)
|
|
fix = anc_score(Fixed(0.5), gen, snr)
|
|
tot["kj"] += kj
|
|
tot["orc"] += orc
|
|
tot["meas"] += meas
|
|
tot["fix"] += fix
|
|
cells += 1
|
|
print(f"{sname:22s}{snr:>6.0f}{kj:>12.2f}{orc:>12.2f}"
|
|
f"{meas:>11.2f}{fix:>11.2f}", flush=True)
|
|
print("-" * 84)
|
|
print(f"{'mean':22s}{'':6s}{tot['kj'] / cells:>12.2f}"
|
|
f"{tot['orc'] / cells:>12.2f}{tot['meas'] / cells:>11.2f}"
|
|
f"{tot['fix'] / cells:>11.2f}")
|
|
survives = tot["meas"] > tot["orc"]
|
|
print(f"\n measured {tot['meas'] / cells:.2f} dB against the oracle "
|
|
f"{tot['orc'] / cells:.2f} dB => "
|
|
f"{'SURVIVES' if survives else 'FAILS'} the objection")
|
|
return survives
|
|
|
|
|
|
if __name__ == "__main__":
|
|
claim_a()
|
|
claim_b()
|