diff --git a/README.md b/README.md index c6fa19b..74d14fd 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,7 @@ platforms, and the same SciPy version scores 2635 for TOMS 748 on Windows agains the same set. The published numbers are from Linux (Ubuntu 24.04, Python 3.12.3). ``` -python -m pytest tests -q # 1240 tests +python -m pytest tests -q # 1205 tests python -u verify_article.py # 548 numbers of the paper against the run data ``` @@ -59,7 +59,7 @@ verify_article.py every number stated in the paper, against the data article/ the paper PREREG_*.md POSTHOC_*.md what was recorded before each run, and what after FINDINGS_*.md what the diagnoses found -tests/ 1240 tests; every bug found in an audit is pinned here +tests/ 1205 tests; every bug found in an audit is pinned here ``` ## What was refuted diff --git a/anc/__init__.py b/anc/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/anc/bench.py b/anc/bench.py deleted file mode 100644 index 22a3723..0000000 --- a/anc/bench.py +++ /dev/null @@ -1,190 +0,0 @@ -"""Active noise control: the step size of an adaptive filter. - -The third domain, and the one where the object we have is most exactly the -object the field argues about. - -In NLMS the coefficient update is - - w <- w + mu * e * x / (||x||^2 + eps) - -and mu = 1 is not an arbitrary scale: it is the step that drives the a -posteriori error at the current sample exactly to zero. Every practical -system uses mu < 1, because zeroing the current sample fits the measurement -noise as well as the signal. So mu is precisely a SAFETY FACTOR on an exact -step - the same object as the safety factor in ODE step control, and the same -object as the projection budget in the root-finding safeguard. - -The whole variable-step-size LMS literature exists because mu is hard to -choose, and it answers with more constants: Kwong and Johnston's rule has -four. What this file asks is whether the quantity mu is guessing can be -measured instead. - - predicted a posteriori error (1 - mu) * e if the model were exact - observed a posteriori error d - w_new . x after the update - -Their ratio says how much of e was signal the filter could actually remove -and how much was noise it could not. That is the safety factor, measured. - -The acoustic paths and the noise regimes are taken verbatim from our earlier -ANC benchmark (branch feature/sem-anc-benchmark of the SEM repository), so -the environment is the established one and not built to suit the result. -""" -from __future__ import annotations - -import math - -import numpy as np - -FS = 8000 -PRIMARY_PATH_LEN = 64 -SECONDARY_PATH_LEN = 32 -FILTER_LEN = 64 - - -# ---------------------------------------------------------------- paths --- -def make_primary_path(rng): - h = np.zeros(PRIMARY_PATH_LEN) - delay = rng.integers(5, 10) - h[delay] = 0.8 + 0.2 * rng.standard_normal() - for k in range(3): - idx = delay + rng.integers(3, 15) - if idx < PRIMARY_PATH_LEN: - h[idx] = (0.3 - 0.1 * k) * rng.standard_normal() - for i in range(delay + 5, PRIMARY_PATH_LEN): - h[i] += 0.05 * math.exp(-0.1 * (i - delay)) * rng.standard_normal() - return h - - -def make_secondary_path(rng): - h = np.zeros(SECONDARY_PATH_LEN) - delay = rng.integers(2, 5) - h[delay] = 0.9 + 0.1 * rng.standard_normal() - for k in range(2): - idx = delay + rng.integers(2, 8) - if idx < SECONDARY_PATH_LEN: - h[idx] = (0.2 - 0.05 * k) * rng.standard_normal() - for i in range(delay + 3, SECONDARY_PATH_LEN): - h[i] += 0.03 * math.exp(-0.15 * (i - delay)) * rng.standard_normal() - return h - - -# ------------------------------------------------------------- regimes --- -def gen_tonal(n, rng, f0=200.0): - t = np.arange(n) / FS - return np.sin(2 * math.pi * f0 * t + rng.uniform(0, 2 * math.pi)) - - -def gen_broadband(n, rng): - pink = np.zeros(n) - rows = np.zeros(5) - for i in range(n): - for j in range(5): - if i % (2 ** j) == 0: - rows[j] = rng.standard_normal() - pink[i] = np.sum(rows) - return pink / (np.std(pink) + 1e-12) - - -def gen_multitonal(n, rng, f0=100.0): - t = np.arange(n) / FS - ph = rng.uniform(0, 2 * math.pi, 3) - sig = (np.sin(2 * math.pi * f0 * t + ph[0]) - + 0.5 * np.sin(2 * math.pi * 2 * f0 * t + ph[1]) - + 0.25 * np.sin(2 * math.pi * 3 * f0 * t + ph[2])) - return sig / (np.max(np.abs(sig)) + 1e-12) - - -def gen_silence(n, rng): - return rng.standard_normal(n) * 0.01 - - -def gen_impulse_train(n, rng): - sig = np.zeros(n) - for _ in range(rng.integers(5, 12)): - pos = rng.integers(0, n) - width = rng.integers(10, 40) - amp = rng.uniform(1.0, 3.0) * rng.choice([-1, 1]) - end = min(pos + width, n) - sig[pos:end] = amp * np.exp(-5.0 * np.arange(end - pos) / width) - return sig - - -def scen_stable(n, rng): - seg = n // 4 - return np.concatenate([gen_tonal(seg, rng), gen_broadband(seg, rng), - gen_multitonal(seg, rng), - gen_silence(n - 3 * seg, rng)]) - - -def scen_gradual(n, rng): - a = np.linspace(0, 1, n) - return (1 - a) * gen_tonal(n, rng) + a * gen_broadband(n, rng) - - -def scen_abrupt(n, rng): - out = np.zeros(n) - gens = [gen_tonal, gen_broadband, gen_multitonal] - pos = 0 - while pos < n: - seg = min(int(rng.integers(n // 20, n // 8)), n - pos) - out[pos:pos + seg] = gens[rng.integers(0, 3)](seg, rng) - pos += seg - return out - - -def scen_cabin(n, rng): - engine = gen_multitonal(n, rng, f0=80.0) - road = 0.6 * gen_broadband(n, rng) - imp = 0.4 * gen_impulse_train(n, rng) - s = engine + road + imp - return s / (np.max(np.abs(s)) + 1e-12) - - -def scen_impulsive(n, rng): - return gen_impulse_train(n, rng) + 0.2 * gen_broadband(n, rng) - - -SCENARIOS = [("stable regimes", scen_stable), ("gradual transition", scen_gradual), - ("abrupt switching", scen_abrupt), ("car cabin", scen_cabin), - ("impulsive", scen_impulsive)] - - -# --------------------------------------------------------- the filter --- -def run_nlms(noise, prim, sec, rule, snr_db=30.0, seed=0): - """One NLMS/FxLMS run. Everything is fixed except `rule`, which returns - the step size for the current sample.""" - rng = np.random.default_rng(seed) - n = len(noise) - d = np.convolve(noise, prim)[:n] # what reaches the ear - xf = np.convolve(noise, sec)[:n] # filtered reference - p = float(np.mean(d * d)) - d = d + rng.standard_normal(n) * math.sqrt(p / (10 ** (snr_db / 10.0))) - - w = np.zeros(FILTER_LEN) - buf = np.zeros(FILTER_LEN) - err = np.zeros(n) - rule.reset() - for i in range(n): - buf[1:] = buf[:-1] - buf[0] = xf[i] - y = float(np.dot(w, buf)) - e = d[i] - y - err[i] = e - nrm = float(np.dot(buf, buf)) + 1e-6 - mu = rule.step(e, buf, nrm) - if mu != 0.0: - w += mu * e * buf / nrm - # the a posteriori error, which is what the rule needs to learn - # how much of e it was able to remove - rule.observe(e, d[i] - float(np.dot(w, buf)), mu) - else: - rule.observe(e, e, 0.0) - # noise reduction: power of the residual against the power of what would - # have reached the ear untreated, measured on the second half so the - # transient is not counted - half = n // 2 - p_in = float(np.mean(d[half:] ** 2)) - p_out = float(np.mean(err[half:] ** 2)) - if p_out <= 0.0 or not math.isfinite(p_out): - return float("nan") - return 10.0 * math.log10(p_in / p_out) diff --git a/anc/rules.py b/anc/rules.py deleted file mode 100644 index 67f66d6..0000000 --- a/anc/rules.py +++ /dev/null @@ -1,204 +0,0 @@ -"""Step-size rules for NLMS. Everything else in the filter is shared. - - mu = 1 drives the a posteriori error of the current sample to zero - mu < 1 hedges, because that sample carries measurement noise too - -So mu is a safety factor on an exact step, exactly as in ODE step control. -The variable-step-size literature answers the difficulty of choosing it with -more constants; the rule at the bottom measures it instead. -""" -from __future__ import annotations - -import math - - -class Rule: - name = "?" - constants = 0 - - def reset(self): - pass - - def step(self, e, x, nrm): - raise NotImplementedError - - def observe(self, e_pre, e_post, mu): - pass - - -class Fixed(Rule): - def __init__(self, mu): - self.mu = mu - self.name = f"fixed mu={mu:g}" - self.constants = 1 - - def step(self, e, x, nrm): - return self.mu - - -class KwongJohnston(Rule): - """The classical variable step size: mu(n+1) = a*mu(n) + g*e(n)^2, clamped. - - Four constants, and the paper that introduced it spends its length on how - to pick them.""" - - name = "Kwong-Johnston VSS" - constants = 4 - A, G, MU_MIN, MU_MAX = 0.97, 0.02, 0.001, 1.0 - - def reset(self): - self.mu = 0.5 - - def step(self, e, x, nrm): - return self.mu - - def observe(self, e_pre, e_post, mu): - self.mu = self.A * self.mu + self.G * e_pre * e_pre - self.mu = min(self.MU_MAX, max(self.MU_MIN, self.mu)) - - -class InSampleMeasured(Rule): - """REFUTED - kept so the mistake is reproducible. - - The first attempt measured the A POSTERIORI error of the sample the update - had just been fitted to. That observable is uninformative by - construction: NLMS with mu = 1 drives it to exactly zero, so the rule - reads "the model was perfect", proposes mu = 1 again, and is locked there. - Measured, it scores 18.50 dB against fixed mu = 1's 18.49 - the same rule - wearing a different name. - - It is the same error as the rho detector in the root-finding work: - measuring a quantity in the one place where it is identically zero.""" - - name = "in-sample (refuted)" - constants = 0 - - def reset(self): - self.rho = 1.0 - - def step(self, e, x, nrm): - return min(1.0, 1.0 / self.rho) if self.rho > 0.0 else 1.0 - - def observe(self, e_pre, e_post, mu): - if mu <= 0.0 or e_pre == 0.0: - return - predicted = abs((1.0 - mu) * e_pre) - observed = abs(e_post) - if predicted <= 0.0: - self.rho = 1.0 + observed / (abs(e_pre) + 1e-300) - else: - r = observed / predicted - if math.isfinite(r) and r > 0.0: - self.rho = math.sqrt(self.rho * max(r, 1e-6)) - - -class Measured(Rule): - """The safety factor measured out of sample. - - What mu is hedging against is the part of the error the filter cannot - cancel - the measurement noise. NLMS theory gives the step that minimises - the mean square deviation as the SIGNAL FRACTION of the error, - - mu = 1 - p_noise / p_error, - - and both quantities are observable. p_error is the running power of e. - p_noise is what e cannot be driven below, so it is estimated as the - smallest error power the run has sustained - the floor the filter keeps - hitting. - - Nothing here is chosen. The two smoothing weights are 1/2, an average of - the current observation with the accumulated one, not a tuned forgetting - factor; and the first samples, having no floor to compare against, take - the exact step - the licence the root-finding safeguard grants. - """ - - name = "measured (no constants)" - constants = 0 - - def reset(self): - self.p_err = None # running error power - self.floor = None # the smallest sustained error power seen - - def step(self, e, x, nrm): - if self.p_err is None or self.floor is None or self.p_err <= 0.0: - return 1.0 - frac = 1.0 - self.floor / self.p_err - return min(1.0, max(0.0, frac)) - - def observe(self, e_pre, e_post, mu): - p = e_pre * e_pre - self.p_err = p if self.p_err is None else 0.5 * (self.p_err + p) - if self.floor is None: - self.floor = self.p_err - elif self.p_err < self.floor: - self.floor = self.p_err - else: - # the floor may rise when the environment changes; let it drift up - # by the same averaging, so a regime change is not remembered for - # ever as a floor that no longer exists - self.floor = 0.5 * (self.floor + min(self.p_err, self.floor * 2.0)) - - -class CosineSimilarity(Rule): - """mu as a SIMILARITY, which is what this project was supposed to be about. - - The two rules above estimate mu from variances - a noise floor, a model - discrepancy. Neither is a similarity; both are statistics of the error. - The similarity formulation asks a different question: - - how much does the error RESEMBLE something the filter can remove? - - The filter can only remove what is correlated with its input. So the - quantity is the cosine similarity between the error and the filtered - reference, - - mu = || / (||e|| ||x||) in [0, 1], - - which is the cosine of the original paper, applied to the two signals the - filter actually has. It equals 1 when the error is entirely explainable - by the reference - the exact step is then right - and 0 when the error is - orthogonal to it, which is measurement noise the update can only fit. - - The three running means use the weight 1/2, an average of the current - observation with the accumulated one, and the first sample takes the exact - step for want of anything to compare. No constant is chosen. - """ - - name = "cosine similarity (no constants)" - constants = 0 - S = None # None means the identity: mu = |cos| directly - - def reset(self): - self.ex = None # running - self.ee = None # running - self.xx = None # running - - def step(self, e, x, nrm): - # decide from what has been accumulated SO FAR, then fold in the - # current sample - so the step never uses the observation it is about - # to be judged on - if self.ex is None: - mu = 1.0 - else: - den = math.sqrt(self.ee * self.xx) - if den <= 0.0 or not math.isfinite(den): - mu = 1.0 - else: - cos = min(1.0, abs(self.ex) / den) - if self.S is None: - mu = cos - else: - # the general family: mu = S(d) with d = 1 - |cos| the - # dissimilarity between the error and what the filter can - # remove. d = 0 means fully removable, so S(0) = 1 and the - # exact step is taken; the identity above is S(d) = 1 - d. - mu = max(0.0, min(1.0, self.S(1.0 - cos))) - xv = float(x[0]) - self.ex = e * xv if self.ex is None else 0.5 * (self.ex + e * xv) - self.ee = e * e if self.ee is None else 0.5 * (self.ee + e * e) - self.xx = xv * xv if self.xx is None else 0.5 * (self.xx + xv * xv) - return mu - - -ALL = [Fixed(0.1), Fixed(0.5), Fixed(1.0), KwongJohnston(), - InSampleMeasured(), Measured(), CosineSimilarity()] diff --git a/odestep/__init__.py b/odestep/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/odestep/controllers.py b/odestep/controllers.py deleted file mode 100644 index b6074be..0000000 --- a/odestep/controllers.py +++ /dev/null @@ -1,157 +0,0 @@ -"""Step-size controllers. Everything else in the integration is SciPy's -(the RK45 tableau, the stage combination, the embedded error estimate), so -these are the only thing under comparison. - -The classical controllers are made of constants somebody chose - a safety -factor, growth clamps, PI gains. The question this module asks is what -happens if the quantity a constant is GUESSING is measured instead. - - root finding step-size control - ------------ ----------------- - the scale of the value axis the safety factor - chosen -> brittle, linear chosen -> works, three constants - measured-> superlinear measured-> works, none - -What did NOT transfer is the credit ledger itself. A first version banked an -accuracy margin exactly as the root-finding safeguard does; the margin turned -out to take no part in any decision and was removed. The move that transfers -is the more general one - replace a chosen constant by the measured value of -what it approximates - not the particular bookkeeping that carried it there. -""" -from __future__ import annotations - -import math - - -class Controller: - name = "?" - constants = 0 - - def reset(self): - pass - - def update(self, h, e, ok, order): - raise NotImplementedError - - -class Elementary(Controller): - """h_new = h * (1/e)^(1/(p+1)) with no safety factor and no clamps. - - Included because it is the controller the error model alone dictates. It - is what every practical code refuses to use, and the reason for that - refusal is exactly what the constants in the others are buying.""" - - name = "elementary (no constants)" - constants = 0 - - def update(self, h, e, ok, order): - if e == 0.0: - return h * 10.0 - return h * (1.0 / e) ** (1.0 / (order + 1)) - - -class Classical(Controller): - """The textbook controller: a safety factor and two growth clamps.""" - - name = "classical (scipy: safety 0.9, clamps 0.2/10)" - constants = 3 - # exactly the constants scipy.integrate._ivp.rk uses, so the driver can - # be validated against solve_ivp with the same control law - SAFETY, MINF, MAXF = 0.9, 0.2, 10.0 - - def update(self, h, e, ok, order): - if e == 0.0: - fac = self.MAXF - else: - fac = self.SAFETY * (1.0 / e) ** (1.0 / (order + 1)) - fac = min(self.MAXF, max(self.MINF, fac)) - return h * fac - - -class GustafssonPI(Controller): - """The PI controller used in modern codes: two gains, a safety factor and - two clamps.""" - - name = "PI (Gustafsson, 5 constants)" - constants = 5 - SAFETY, MINF, MAXF = 0.9, 0.2, 10.0 - K_I, K_P = 0.3, 0.4 - - def reset(self): - self.e_prev = None - - def update(self, h, e, ok, order): - if e == 0.0: - return h * self.MAXF - p1 = order + 1 - if self.e_prev is None or not ok: - fac = self.SAFETY * (1.0 / e) ** (1.0 / p1) - else: - fac = (self.SAFETY * (1.0 / e) ** (self.K_I / p1) - * (self.e_prev / e) ** (self.K_P / p1)) - fac = min(self.MAXF, max(self.MINF, fac)) - if ok: - self.e_prev = e - return h * fac - - -class MeasuredSafety(Controller): - """The safety factor measured instead of chosen. - - NOTE ON THE NAME. This began as "the credit rule transferred", and it - banked a margin exactly as the root-finding safeguard does. The margin - turned out to play no part in the decision - the step factor was computed - from the measured model discrepancy alone - so the banking was vestigial - and the name was unearned. What transfers is not the credit ledger but - the deeper move behind it: replace a chosen constant with the quantity it - is guessing, measured from the run. The ledger is left out because it did - nothing here, and calling this a credit controller would be a mislabel. - - The elementary controller aims at err = tol exactly, and is refused in - practice because the model err ~ C h^(p+1) is only approximate: aiming at - the boundary means missing it about half the time, and every miss costs a - whole rejected step. The classical answer is a safety factor chosen by - hand. - - The answer here is to measure the thing the safety factor is guessing. - After proposing a step, the model predicts what its error will be; the - step is taken and the error is observed. The discrepancy - - delta = log2(observed / predicted) - - is exactly how wrong the model was, in bits, and the next step aims at - 2^-delta instead of at 1. When the model is right, delta -> 0 and the - controller becomes the elementary one; when the model has just been - optimistic, the controller backs off by precisely the observed amount. - - The first step has nothing to have learned from, so it aims at the - tolerance itself - the same licence the root-finding guard grants. - """ - - name = "measured safety (no constants)" - constants = 0 - - def reset(self): - self.predicted = None # error the model expected for this step - self.delta = 0.0 # bits by which the model was last wrong - - def update(self, h, e, ok, order): - p1 = order + 1 - if self.predicted is not None and self.predicted > 0.0 and e > 0.0: - # only optimism matters: the model predicting too SMALL an error - # is what causes rejections. Pessimism costs nothing and is not - # corrected for. - self.delta = max(0.0, math.log2(e / self.predicted)) - if e == 0.0: - fac = 10.0 - else: - # aim at 2^-delta of the tolerance rather than at the tolerance - fac = (2.0 ** (-self.delta) / e) ** (1.0 / p1) - h_new = h * fac - # record what the model now predicts for the step just chosen, so the - # next call can measure how wrong it was - self.predicted = e * (h_new / h) ** p1 - return h_new - - -ALL = [Elementary(), Classical(), GustafssonPI(), MeasuredSafety()] diff --git a/odestep/driver.py b/odestep/driver.py deleted file mode 100644 index ef5a8f5..0000000 --- a/odestep/driver.py +++ /dev/null @@ -1,114 +0,0 @@ -"""Integration driver whose only original part is the step-size controller. - -The tableau, the stage combination and the embedded error estimate are -SciPy's - `scipy.integrate._ivp.rk.rk_step` with RK45's A, B, C and E - so -nothing about the method itself is a transcription of ours. An earlier -version of this file carried a hand-typed Dormand-Prince tableau; that was an -unnecessary risk and an unnecessary thing for a reader to have to check. - -What remains ours is the loop that decides the next step size, which is the -object under study. Everything else is shared between controllers, so a -difference between two runs can only come from the controller. - -The driver is validated against SciPy's own `solve_ivp(method="RK45")`: with -the classical controller - SciPy's own safety factor and clamps - the two -must agree closely on the number of right-hand-side evaluations. -""" -from __future__ import annotations - -import math - -import numpy as np -from scipy.integrate._ivp.rk import RK45, rk_step - -ORDER = RK45.error_estimator_order # 4; the step exponent is ORDER + 1 - - -class Counter: - """Counts evaluations of the right-hand side - the cost model.""" - - __slots__ = ("f", "n") - - def __init__(self, f): - self.f = f - self.n = 0 - - def __call__(self, t, y): - self.n += 1 - return np.asarray(self.f(t, y), dtype=float) - - -def _norm(v, y, y_new, atol, rtol): - scale = atol + rtol * np.maximum(np.abs(y), np.abs(y_new)) - return float(np.sqrt(np.mean((v / scale) ** 2))) - - -def integrate(f, t0, t1, y0, controller, atol=1e-8, rtol=1e-8, - max_steps=200000): - """Integrate with the given controller. Returns what the comparison - needs. The initial step is chosen the same way for every controller, so - none can win by guessing it better.""" - g = Counter(f) - t = float(t0) - y = np.asarray(y0, dtype=float) - fv = g(t, y) - K = np.empty((RK45.n_stages + 1, y.size), dtype=float) - - d0 = _norm(y, y, y, atol, rtol) - d1 = _norm(fv, y, y, atol, rtol) - h = 1e-6 if (d0 < 1e-5 or d1 < 1e-5) else 0.01 * d0 / d1 - h = min(h, abs(t1 - t0)) - controller.reset() - - accepted = rejected = 0 - while t < t1 and accepted + rejected < max_steps: - h = min(h, t1 - t) - if h <= 0.0: - break - y_new, f_new = rk_step(g, t, y, fv, h, RK45.A, RK45.B, RK45.C, K) - err = np.dot(K.T, RK45.E) * h # RK45's own error estimate - e = _norm(err, y, y_new, atol, rtol) - ok = e <= 1.0 - h_next = controller.update(h, e, ok, ORDER) - if ok: - accepted += 1 - t, y, fv = t + h, y_new, f_new - else: - rejected += 1 - if not math.isfinite(h_next) or h_next <= 0.0: - break - h = h_next - return {"t": t, "y": list(y), "nfev": g.n, "accepted": accepted, - "rejected": rejected, - "reached_end": abs(t - t1) < 1e-12 * max(1.0, abs(t1))} - - -def validate_against_scipy(atol=1e-8, rtol=1e-8): - """The driver against SciPy's own solver, with SciPy's own controller. - - This replaces the hand-written order check that the previous version - needed: if the driver reproduces solve_ivp's evaluation count with the - same control law, the tableau, the stage combination, the error estimate - and the loop are all the library's behaviour.""" - from scipy.integrate import solve_ivp - - from .controllers import Classical - - out = [] - probs = [ - ("decay", lambda t, y: [-y[0]], 0.0, 10.0, [1.0]), - ("forced", lambda t, y: [-y[0] + math.sin(t)], 0.0, 10.0, [0.0]), - ("lotka", lambda t, y: [1.5 * y[0] - y[0] * y[1], - -3.0 * y[1] + y[0] * y[1]], 0.0, 15.0, - [10.0, 5.0]), - ("vdp", lambda t, y: [y[1], 5.0 * (1 - y[0] ** 2) * y[1] - y[0]], - 0.0, 25.0, [2.0, 0.0]), - ] - for name, f, t0, t1, y0 in probs: - mine = integrate(f, t0, t1, y0, Classical(), atol=atol, rtol=rtol) - ref = solve_ivp(f, (t0, t1), y0, method="RK45", atol=atol, - rtol=rtol) - out.append((name, mine["nfev"], int(ref.nfev), - max(abs(a - b) for a, b in zip(mine["y"], - ref.y[:, -1])))) - return out diff --git a/relax/__init__.py b/relax/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/relax/sor.py b/relax/sor.py deleted file mode 100644 index 9b828b7..0000000 --- a/relax/sor.py +++ /dev/null @@ -1,238 +0,0 @@ -"""Relaxation for iterative linear solvers, with the relaxation factor -pluggable. - -The relaxation parameter omega of SOR is the oldest hand-tuned constant in -numerical linear algebra. Convergence requires 0 < omega < 2; the optimal -value is known in closed form only for model problems, and everywhere else it -is guessed, swept, or left at a folklore value - the under-relaxation factors -of 0.7 and 0.3 in every SIMPLE-type CFD code being the best-known example. - -The similarity formulation. SOR over-relaxes because the Gauss-Seidel -correction under-shoots CONSISTENTLY IN THE SAME DIRECTION: the iteration -creeps towards the solution instead of stepping to it. How consistent the -direction is, is a similarity between successive corrections: - - s = cos(d_k, d_{k-1}) in [-1, 1], omega = 1 + max(0, s) in [1, 2] - -s = 1 means the correction is repeating itself exactly, so the step should be -doubled - the maximum over-relaxation the theory allows. s <= 0 means the -corrections are fighting each other, so no over-relaxation is warranted and -the method falls back to Gauss-Seidel. The admissible interval of the theory -and the range of the similarity coincide, and nothing is chosen. - -The sweeps below are written out because SciPy ships no SOR. They are -validated against theory: for the model Poisson problem the SOR spectral -radius is exactly omega_opt - 1, which the measured convergence rate must -reproduce. -""" -from __future__ import annotations - -import math - -import numpy as np - - -def poisson_rhs(n, rng=None): - """A right-hand side on an n x n interior grid.""" - if rng is None: - rng = np.random.default_rng(0) - return rng.standard_normal((n, n)) - - -def sweep_gs(u, f, h2, omega, coeff=None): - """One Gauss-Seidel / SOR sweep in place, returning the correction. - - coeff is None for the 5-point Laplacian; otherwise (cw, ce, cs, cn, cc) - arrays giving a variable-coefficient or anisotropic stencil. - """ - n = u.shape[0] - before = u.copy() - if coeff is None: - for i in range(n): - for j in range(n): - s = 0.0 - if i > 0: - s += u[i - 1, j] - if i < n - 1: - s += u[i + 1, j] - if j > 0: - s += u[i, j - 1] - if j < n - 1: - s += u[i, j + 1] - new = (s - h2 * f[i, j]) / 4.0 - u[i, j] += omega * (new - u[i, j]) - else: - cw, ce, cs, cn, cc = coeff - for i in range(n): - for j in range(n): - s = 0.0 - if i > 0: - s += cw[i, j] * u[i - 1, j] - if i < n - 1: - s += ce[i, j] * u[i + 1, j] - if j > 0: - s += cs[i, j] * u[i, j - 1] - if j < n - 1: - s += cn[i, j] * u[i, j + 1] - new = (s - h2 * f[i, j]) / cc[i, j] - u[i, j] += omega * (new - u[i, j]) - return u - before - - -def residual(u, f, h2, coeff=None): - n = u.shape[0] - r = np.zeros_like(u) - if coeff is None: - for i in range(n): - for j in range(n): - s = -4.0 * u[i, j] - if i > 0: - s += u[i - 1, j] - if i < n - 1: - s += u[i + 1, j] - if j > 0: - s += u[i, j - 1] - if j < n - 1: - s += u[i, j + 1] - r[i, j] = h2 * f[i, j] - s - else: - cw, ce, cs, cn, cc = coeff - for i in range(n): - for j in range(n): - s = -cc[i, j] * u[i, j] - if i > 0: - s += cw[i, j] * u[i - 1, j] - if i < n - 1: - s += ce[i, j] * u[i + 1, j] - if j > 0: - s += cs[i, j] * u[i, j - 1] - if j < n - 1: - s += cn[i, j] * u[i, j + 1] - r[i, j] = h2 * f[i, j] - s - return r - - -class Fixed: - def __init__(self, omega): - self.omega = omega - self.name = f"omega = {omega:g}" - self.constants = 1 - - def reset(self): - pass - - def next_omega(self, d_prev, d_now): - return self.omega - - -class OptimalPoisson(Fixed): - """The closed-form optimum for the model problem - an oracle that needs no - tuning because the answer is a theorem.""" - - def __init__(self, n): - w = 2.0 / (1.0 + math.sin(math.pi / (n + 1))) - super().__init__(w) - self.name = f"omega* = {w:.4f} (theory)" - - -class SimilarityOmega: - """omega = 1 + max(0, cos(d_k, d_{k-1})). - - The first sweep has no predecessor to compare with, so it runs at - Gauss-Seidel - the same licence granted in the root-finding safeguard. - """ - - name = "similarity (no constants)" - constants = 0 - - def reset(self): - pass - - def next_omega(self, d_prev, d_now): - if d_prev is None or d_now is None: - return 1.0 - a = float(np.vdot(d_prev, d_now)) - na = float(np.linalg.norm(d_prev)) * float(np.linalg.norm(d_now)) - if na <= 0.0 or not math.isfinite(na): - return 1.0 - return 1.0 + max(0.0, a / na) - - -class ContractionOmega: - """omega driven by the contraction the iteration actually achieves. - - The direction-similarity rule above over-relaxes: the cosine between - successive SOR corrections sits near 1 almost always, so 1 + cos saturates - at 2. Measured, it picks 1.87 where the optimum is 1.777 and 1.43 where - it is 1.2. The measure was right and the map into omega was not. - - The map that is derived rather than chosen comes from the theory itself: - for omega >= omega* the SOR spectral radius is exactly omega - 1. So the - observed contraction of the residual, rho = ||r_k|| / ||r_{k-1}||, says on - which side of the optimum the iteration is sitting: - - rho > omega - 1 -> below the optimum, lengthen - rho ~ omega - 1 -> at or above it, hold - - and the correction is the discrepancy itself, - - omega <- omega + (rho - (omega - 1)), clipped to [1, 2). - - The clip is the convergence condition of the method, not a tuning range, - and the first sweep runs at omega = 1 for want of a previous residual. - - REFUTED, and worse than the rule it was meant to repair. Measured: 235 - sweeps against the direction rule's 107 on the model problem, 589 against - 102 on the anisotropic one. The update is undamped, so on the first - sweep, with omega = 1 and rho near 1, it sets omega = 1 + (1 - 0) = 2 and - stays pinned at the limit. Damping it would take a constant, which is the - thing the exercise is trying to avoid, so the line stops here rather than - trying a third map. - """ - - name = "contraction-driven (no constants)" - constants = 0 - - def reset(self): - self.omega = 1.0 - self.r_prev = None - - def next_omega(self, d_prev, d_now): - return self.omega - - def note_residual(self, rn): - if self.r_prev is not None and self.r_prev > 0.0: - rho = rn / self.r_prev - if math.isfinite(rho): - w = self.omega + (rho - (self.omega - 1.0)) - # 2 is the convergence limit of SOR, 1 is Gauss-Seidel - self.omega = min(1.999, max(1.0, w)) - self.r_prev = rn - - -def solve(f, rule, tol=1e-8, maxit=20000, coeff=None): - """Iterate to ||r|| <= tol * ||r_0||. Cost is sweeps, which is what a - practitioner pays. Everything except `rule` is shared.""" - n = f.shape[0] - h2 = 1.0 / (n + 1) ** 2 - u = np.zeros_like(f) - r0 = float(np.linalg.norm(residual(u, f, h2, coeff))) - if r0 == 0.0: - return 0, 0.0 - rule.reset() - d1 = d2 = None # the two most recent corrections, newest first - omegas = [] - for k in range(1, maxit + 1): - # the rule sees the last two corrections; on the first two sweeps one - # of them does not exist yet and the rule returns 1 (Gauss-Seidel) - omega = rule.next_omega(d2, d1) - omegas.append(omega) - d1, d2 = sweep_gs(u, f, h2, omega, coeff), d1 - rn = float(np.linalg.norm(residual(u, f, h2, coeff))) - if hasattr(rule, "note_residual"): - rule.note_residual(rn) - if not math.isfinite(rn): - return maxit, float("nan") - if rn <= tol * r0: - return k, float(np.mean(omegas)) - return maxit, float(np.mean(omegas)) diff --git a/stage15_ode.py b/stage15_ode.py deleted file mode 100644 index c3f27cc..0000000 --- a/stage15_ode.py +++ /dev/null @@ -1,189 +0,0 @@ -"""Stage 15 - the credit rule in a second domain: ODE step-size control. - -The root-finding work ended with one transferable object, the credit -safeguard, and one honest conclusion: as a one-dimensional root finder the -method is subsumed by Chandrupatla 1997. Adaptive step-size control is the -natural place to test whether the transferable object transfers, because it -has the same shape - propose a bold step, look at the outcome, adjust - and -because every practical controller in it is made of constants somebody chose. - - classical scipy's own: safety 0.9, clamps 0.2 and 10 3 constants - PI the same plus two gains 5 constants - elementary what the error model alone dictates 0 constants - measured the safety factor MEASURED instead of chosen 0 constants - -Everything except the controller comes from SciPy: the RK45 tableau, the -stage combination and the embedded error estimate, via -scipy.integrate._ivp.rk. The driver is checked against solve_ivp(RK45) with -the same control law before any comparison is made. A difference between two -rows can only come from the controller. - -Cost model: evaluations of the right-hand side, as evaluations of f were in -the root-finding work. A run that does not reach the end of the interval is -a failure and is reported as one. Accuracy is checked against a reference so -that a controller cannot win by being fast and wrong. -""" -from __future__ import annotations - -import math - -from odestep.controllers import ALL -from odestep.driver import integrate, validate_against_scipy - -TOLS = (1e-4, 1e-6, 1e-8, 1e-10) - - -# -------------------------------------------------------------------------- -# problems: classics with known or high-accuracy reference endpoints -# -------------------------------------------------------------------------- -def p_decay(): - f = lambda t, y: [-y[0]] - return "exponential decay", f, 0.0, 10.0, [1.0], [math.exp(-10.0)] - - -def p_forced(): - f = lambda t, y: [-y[0] + math.sin(t)] - ex = 0.5 * (math.sin(10.0) - math.cos(10.0) + math.exp(-10.0)) - return "y' = -y + sin t", f, 0.0, 10.0, [0.0], [ex] - - -def p_lotka(): - def f(t, y): - return [1.5 * y[0] - y[0] * y[1], -3.0 * y[1] + y[0] * y[1]] - return "Lotka-Volterra", f, 0.0, 15.0, [10.0, 5.0], None - - -def p_vanderpol(mu=5.0): - def f(t, y): - return [y[1], mu * (1 - y[0] * y[0]) * y[1] - y[0]] - return f"van der Pol mu={mu:g}", f, 0.0, 25.0, [2.0, 0.0], None - - -def p_brusselator(): - def f(t, y): - return [1.0 + y[0] * y[0] * y[1] - 4.0 * y[0], - 3.0 * y[0] - y[0] * y[0] * y[1]] - return "Brusselator", f, 0.0, 20.0, [1.5, 3.0], None - - -def p_kepler(ecc): - """The two-body problem at eccentricity e - DETEST's D-class, the classic - hard non-stiff test: the solution has a sharp passage at perihelion.""" - def f(t, y): - x1, x2, v1, v2 = y - r3 = (x1 * x1 + x2 * x2) ** 1.5 - return [v1, v2, -x1 / r3, -x2 / r3] - y0 = [1.0 - ecc, 0.0, 0.0, math.sqrt((1 + ecc) / (1 - ecc))] - return f"Kepler e={ecc:g}", f, 0.0, 20.0, y0, None - - -def p_arenstorf(): - """The restricted three-body orbit, the standard showcase for adaptive - step control: the step must vary by orders of magnitude along the path.""" - mu = 0.012277471 - mp = 1.0 - mu - - def f(t, y): - y1, y2, y3, y4 = y - d1 = ((y1 + mu) ** 2 + y2 * y2) ** 1.5 - d2 = ((y1 - mp) ** 2 + y2 * y2) ** 1.5 - return [y3, y4, - y1 + 2 * y4 - mp * (y1 + mu) / d1 - mu * (y1 - mp) / d2, - y2 - 2 * y3 - mp * y2 / d1 - mu * y2 / d2] - y0 = [0.994, 0.0, 0.0, -2.00158510637908252240537862224] - return "Arenstorf orbit", f, 0.0, 17.0652165601579625588917206249, y0, y0 - - -PROBLEMS = [p_decay(), p_forced(), p_lotka(), p_vanderpol(), p_brusselator(), - p_kepler(0.5), p_kepler(0.9), p_arenstorf()] - - -def reference(f, t0, t1, y0): - """A reference endpoint from the same integrator at a far tighter - tolerance - not from a different code, so that a disagreement cannot be - blamed on a different method.""" - from odestep.controllers import Classical - r = integrate(f, t0, t1, y0, Classical(), atol=1e-13, rtol=1e-13, - max_steps=2000000) - return r["y"] if r["reached_end"] else None - - -def scaled_err(y, ref, atol, rtol): - """The achieved error in the SAME mixed norm the integrator controls. - - A purely relative criterion is wrong here and was: the decay problem ends - at exp(-10) ~ 4.5e-5, so at atol = 1e-4 the tolerance permits a hundred - per cent relative error, and a correct run was being marked as a failure. - What a run must be judged against is the accuracy it was asked for.""" - if ref is None: - return float("nan") - s = 0.0 - for a, b in zip(y, ref): - sc = atol + rtol * abs(b) - s += ((a - b) / sc) ** 2 - return math.sqrt(s / len(ref)) - - -def main(): - print("instrument check: the driver against scipy solve_ivp(RK45), " - "same control law") - worst = 0.0 - for name, mine, ref, d in validate_against_scipy(): - rel = abs(mine - ref) / ref - worst = max(worst, rel) - print(f" {name:10s} ours {mine:>6d} solve_ivp {ref:>6d} " - f"{100 * (mine - ref) / ref:+6.1f}% endpoint agrees to {d:.1e}") - if worst > 0.05: - print(" driver does NOT reproduce the library - stopping") - return - print(" the tableau, the stages and the error estimate are scipy's;") - print(" only the controller is ours\n") - - refs = {} - for name, f, t0, t1, y0, exact in PROBLEMS: - refs[name] = exact if exact is not None else reference(f, t0, t1, y0) - - grand = {c.name: [0, 0, 0] for c in ALL} # nfev, rejected, failures - for tol in TOLS: - print("=" * 96) - print(f"tolerance {tol:g} (atol = rtol = tol)") - print("=" * 96) - print(f"{'problem':22s}" + "".join(f"{c.name.split(' ')[0]:>15s}" - for c in ALL)) - print("-" * 96) - for name, f, t0, t1, y0, _exact in PROBLEMS: - cells = [] - for c in ALL: - r = integrate(f, t0, t1, y0, c, atol=tol, rtol=tol) - err = scaled_err(r["y"], refs[name], tol, tol) - # local error control does not bound the global error, so the - # bar is generous: a run fails only if it misses the accuracy - # it asked for by more than three orders of magnitude, or does - # not reach the end at all - bad = (not r["reached_end"]) or ( - math.isfinite(err) and err > 1000.0) - grand[c.name][0] += r["nfev"] - grand[c.name][1] += r["rejected"] - grand[c.name][2] += bad - cells.append(f"{r['nfev']:>9d}{'!' if bad else ' '}" - f"{r['rejected']:>4d}") - print(f"{name:22s}" + "".join(f"{c:>15s}" for c in cells)) - print() - - print("=" * 96) - print("TOTALS over every problem and tolerance " - "(nfev, rejected steps, failures)") - print("=" * 96) - print(f"{'controller':34s}{'constants':>11s}{'nfev':>10s}" - f"{'rejected':>11s}{'failures':>11s}") - print("-" * 96) - for c in ALL: - n, rej, bad = grand[c.name] - print(f"{c.name:34s}{c.constants:>11d}{n:>10d}{rej:>11d}{bad:>11d}") - print("\n'!' marks a run that did not reach the end or missed the " - "reference by more\nthan a thousand tolerances - a controller " - "cannot win by being fast and wrong.") - - -if __name__ == "__main__": - main() diff --git a/stage17_anc.py b/stage17_anc.py deleted file mode 100644 index 39e5468..0000000 --- a/stage17_anc.py +++ /dev/null @@ -1,81 +0,0 @@ -"""Stage 17 - the third domain: the step size of an adaptive noise canceller. - -Root finding: the method lands second behind Chandrupatla 1997. -ODE step control: level with a five-constant PI controller, behind the -three-constant classical one. -Softening somebody else's acceptance test: fails, three times over. - -Active noise control is the domain where the object we have - a step size -that is a safety factor on an exact step - is exactly the object the field -argues about, and where the standard answers are openly constant-ridden. - -The acoustic paths, the noise regimes and the scenarios come from our earlier -ANC benchmark (SEM repository, branch feature/sem-anc-benchmark), so the -environment is one that existed before this question was asked. - -Metric: noise reduction in dB on the second half of each run, so that the -adaptation transient is not counted. Higher is better. Ten seeds per cell. -""" -from __future__ import annotations - -import math - -import numpy as np - -from anc.bench import (FS, SCENARIOS, make_primary_path, make_secondary_path, - run_nlms) -from anc.rules import ALL - -N = 8000 # one second -SEEDS = 10 -SNRS = (30.0, 15.0) - - -def main(): - print(f"{len(SCENARIOS)} scenarios x {len(SNRS)} noise levels x {SEEDS} " - f"seeds, {N / FS:g} s each\n") - grand = {r.name: [] for r in ALL} - for snr in SNRS: - print("=" * 88) - print(f"reference SNR {snr:g} dB noise reduction, dB (higher is " - f"better)") - print("=" * 88) - print(f"{'scenario':22s}" + "".join(f"{r.name[:13]:>13s}" for r in ALL)) - print("-" * 88) - for sname, gen in SCENARIOS: - row = {r.name: [] for r in ALL} - for s in range(SEEDS): - rng = np.random.default_rng(1000 + s) - noise = gen(N, rng) - prim = make_primary_path(rng) - sec = make_secondary_path(rng) - for r in ALL: - nr = run_nlms(noise, prim, sec, r, snr_db=snr, seed=s) - if math.isfinite(nr): - row[r.name].append(nr) - grand[r.name].append(nr) - cells = [] - for r in ALL: - v = row[r.name] - cells.append(f"{(sum(v) / len(v)) if v else float('nan'):>13.2f}") - print(f"{sname:22s}" + "".join(cells), flush=True) - print() - - print("=" * 88) - print("OVERALL") - print("=" * 88) - print(f"{'rule':26s}{'constants':>11s}{'mean NR dB':>13s}" - f"{'worst cell':>13s}{'diverged':>11s}") - print("-" * 88) - for r in ALL: - v = grand[r.name] - n_bad = SEEDS * len(SCENARIOS) * len(SNRS) - len(v) - print(f"{r.name:26s}{r.constants:>11d}" - f"{(sum(v) / len(v)) if v else float('nan'):>13.2f}" - f"{min(v) if v else float('nan'):>13.2f}{n_bad:>11d}") - print("\n'diverged' counts runs whose residual power was not finite - the " - "failure mode\na step size that is too large actually produces.") - - -if __name__ == "__main__": - main() diff --git a/stage18_bulletproof.py b/stage18_bulletproof.py index f93d281..35ac037 100644 --- a/stage18_bulletproof.py +++ b/stage18_bulletproof.py @@ -12,15 +12,6 @@ falsi goes from 133 solved and 18799 evaluations to 193 and 8056. 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 @@ -32,8 +23,6 @@ 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 @@ -162,69 +151,5 @@ def claim_a(): 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() diff --git a/stage20_prediction.py b/stage20_prediction.py deleted file mode 100644 index c1e4f49..0000000 --- a/stage20_prediction.py +++ /dev/null @@ -1,204 +0,0 @@ -"""Stage 20 - a prediction registered before the run, then the run. - -The request was to find the domain where the method wins unambiguously. The -way to get that wrong is to try domains until one of them comes out ahead; -that is fitting, spread over time. So this file states the prediction and the -falsification condition FIRST, from what the quantity does, and then measures. - -WHAT MU IS. In NLMS the update w <- w + mu * e * x / ||x||^2 with mu = 1 -drives the a posteriori error of the current sample to exactly zero. Every -practical system uses mu < 1 because part of e is measurement noise the filter -cannot cancel, and fitting it is harmful. So mu is the SIGNAL FRACTION of the -error, and its correct value is - - mu* = 1 - p_noise / p_error. - -THE PREDICTION. mu* is a property of the noise level. Therefore: - - (A) CONTROL - when the noise level is constant, a constant mu can be exactly - right. An oracle that picks the best fixed mu for that scenario should - MATCH OR BEAT the measured rule. If the measured rule wins here, the - reasoning is wrong and everything below is suspect. - - (B) PREDICTION - when the noise level CHANGES during the run, no constant is - right throughout, while a measured quantity tracks. The measured rule - should BEAT the oracle fixed mu, and the margin should grow with the size - of the change. - -FALSIFICATION. If (B) fails - if an oracle constant matches the measured rule -even when the noise level swings by 25 dB within a single run - then there is -no such domain on this evidence, and that is the finding. - -The oracle is the strongest competitor available: the best of a mu grid chosen -separately for every scenario, with knowledge of that scenario, which no -deployable controller has. Kwong-Johnston is given the same treatment. -""" -from __future__ import annotations - -import math - -import numpy as np - -from anc.bench import (FS, gen_broadband, gen_impulse_train, gen_multitonal, - gen_tonal, make_primary_path, make_secondary_path) -from anc.rules import Fixed, KwongJohnston, Measured, Rule - -N = 16000 # two seconds -SEEDS = 8 -MU_GRID = (0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.85, 1.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)] -FILTER_LEN = 64 - - -class KJTuned(KwongJohnston): - def __init__(self, a, g): - self.A, self.G = a, g - self.name = f"KJ({a:g},{g:g})" - - -def run(noise, prim, sec, rule, snr_profile, seed=0): - """One run. snr_profile(i) gives the instantaneous SNR in dB, so the - stationary and the varying case go through exactly the same code.""" - rng = np.random.default_rng(seed) - n = len(noise) - d = np.convolve(noise, prim)[:n] - xf = np.convolve(noise, sec)[:n] - - # the disturbance is added with a time-varying power, which is the only - # thing that differs between the control and the prediction - p_sig = float(np.mean(d * d)) + 1e-30 - snr = np.array([snr_profile(i / n) for i in range(n)]) - amp = np.sqrt(p_sig / (10.0 ** (snr / 10.0))) - d = d + rng.standard_normal(n) * amp - - w = np.zeros(FILTER_LEN) - buf = np.zeros(FILTER_LEN) - err = np.zeros(n) - rule.reset() - for i in range(n): - buf[1:] = buf[:-1] - buf[0] = xf[i] - e = d[i] - float(np.dot(w, buf)) - err[i] = e - nrm = float(np.dot(buf, buf)) + 1e-6 - mu = rule.step(e, buf, nrm) - if mu != 0.0: - w += mu * e * buf / nrm - rule.observe(e, d[i] - float(np.dot(w, buf)), mu) - else: - rule.observe(e, e, 0.0) - # The measurement window matters, and the first version got it wrong. - # Scoring only the second half made the varying-noise scenarios - # indistinguishable from a constant run at their FINAL level - the ramp - # and the step have both arrived by then, so the metric never saw the - # transition and the prediction could not have been confirmed even if it - # were true. The window now excludes only the initial adaptation - # transient, which every controller pays alike. - warm = n // 10 - p_in = float(np.mean(d[warm:] ** 2)) - p_out = float(np.mean(err[warm:] ** 2)) - if p_out <= 0.0 or not math.isfinite(p_out): - return float("nan") - return 10.0 * math.log10(p_in / p_out) - - -def score(rule, gen, profile, seeds=SEEDS): - vals = [] - for s in range(seeds): - rng = np.random.default_rng(2000 + s) - noise = gen(N, rng) - nr = run(noise, make_primary_path(rng), make_secondary_path(rng), - rule, profile, seed=s) - if math.isfinite(nr): - vals.append(nr) - return (sum(vals) / len(vals)) if vals else float("nan") - - -def oracle_fixed(gen, profile): - """The ceiling of what any constant can do on this scenario: the best mu - of the grid, chosen WITH knowledge of the scenario.""" - best, arg = -1e30, None - for mu in MU_GRID: - v = score(Fixed(mu), gen, profile) - if math.isfinite(v) and v > best: - best, arg = v, mu - return best, arg - - -def oracle_kj(gen, profile): - best, arg = -1e30, None - for a, g in KJ_GRID: - v = score(KJTuned(a, g), gen, profile) - if math.isfinite(v) and v > best: - best, arg = v, (a, g) - return best, arg - - -# --------------------------------------------------------------- scenarios -- -def sig_tonal(n, rng): - return gen_tonal(n, rng) - - -def sig_cabin(n, rng): - s = gen_multitonal(n, rng, f0=80.0) + 0.6 * gen_broadband(n, rng) \ - + 0.4 * gen_impulse_train(n, rng) - return s / (np.max(np.abs(s)) + 1e-12) - - -SIGNALS = [("tonal", sig_tonal), ("car cabin", sig_cabin)] - -# noise-level profiles, as functions of the fraction of the run elapsed -FLAT_30 = lambda u: 30.0 -FLAT_15 = lambda u: 15.0 -FLAT_5 = lambda u: 5.0 -SWING_10 = lambda u: 25.0 - 10.0 * u -SWING_25 = lambda u: 30.0 - 25.0 * u -STEP_25 = lambda u: 30.0 if u < 0.5 else 5.0 - -PROFILES = [("flat 30 dB", FLAT_30, "control"), - ("flat 15 dB", FLAT_15, "control"), - ("flat 5 dB", FLAT_5, "control"), - ("ramp 25->15 dB", SWING_10, "prediction"), - ("ramp 30->5 dB", SWING_25, "prediction"), - ("step 30->5 dB", STEP_25, "prediction")] - - -def main(): - print(__doc__.split("The oracle is")[0]) - print("=" * 92) - print("THE RUN") - print("=" * 92) - print(f"{'signal':11s}{'noise profile':17s}{'kind':12s}" - f"{'oracle mu':>11s}{'best mu':>9s}{'oracle KJ':>11s}" - f"{'measured':>10s}{'verdict':>12s}") - print("-" * 92) - control_wins = pred_wins = control_n = pred_n = 0 - for signame, gen in SIGNALS: - for pname, prof, kind in PROFILES: - om, arg = oracle_fixed(gen, prof) - okj, _a = oracle_kj(gen, prof) - ms = score(Measured(), gen, prof) - beats = ms > max(om, okj) - if kind == "control": - control_n += 1 - control_wins += beats - else: - pred_n += 1 - pred_wins += beats - print(f"{signame:11s}{pname:17s}{kind:12s}{om:>11.2f}" - f"{arg:>9.2f}{okj:>11.2f}{ms:>10.2f}" - f"{'MEASURED' if beats else 'oracle':>12s}", flush=True) - print("-" * 92) - print(f"\n control cells (constant noise): measured wins " - f"{control_wins} of {control_n}") - print(f" predicted cells (varying noise): measured wins " - f"{pred_wins} of {pred_n}") - ok = control_wins == 0 and pred_wins == pred_n - print(f"\n the prediction is {'CONFIRMED' if ok else 'NOT confirmed'}: " - f"it required the measured rule to lose every\n control cell and " - f"win every predicted one.") - - -if __name__ == "__main__": - main() diff --git a/stage21_cosine_vs_oracle.py b/stage21_cosine_vs_oracle.py deleted file mode 100644 index 4961173..0000000 --- a/stage21_cosine_vs_oracle.py +++ /dev/null @@ -1,131 +0,0 @@ -"""Stage 21 - the WHOLE similarity family for mu, against the oracle constant. - -Why this is not another throw of the dice. The first two constant-free rules -estimated mu from VARIANCES - a noise floor, a model discrepancy. Neither is -a similarity, and this project is about similarity. The formulation the -thesis dictates asks how much the error RESEMBLES what the filter can remove, -and the filter can only remove what correlates with its input. With - - d = 1 - || / (||e|| ||x||) the dissimilarity, 0 = removable - -the family is mu = S(d) over the same nine similarity functions the -root-finding work used. Testing one member would be as incomplete here as it -would have been there. - -A PREDICTION, FROM THE LAW OF THE FIRST DOMAIN. The contraction law found in -root finding says the shape of S is nearly irrelevant and the slope of (1 - S) -at the origin decides everything: functions approaching 1 linearly work, -functions approaching it quadratically do not. If that law is a property of -similarity rules and not of root finding, it must order this family too: - - predicted strong 1 - d, e^-d, 1/(1+d) slope 1 - predicted middle (1+d)^-2, 1 - sqrt(d), 1 - d^(1/3) - predicted weak 1/(1+d^2), e^-(d^2), 1/sqrt(1+d^2) slope 0 - -THE BAR FOR THE METHOD, STATED BEFORE THE RUN. The competitor is the -strongest a constant can be: the best mu of a grid chosen with knowledge of -the scenario and the noise profile, which no deployable controller has. For -the claim to be bullet-proof the best similarity must beat that oracle on a -MAJORITY of the twelve cells, stationary cells included, where a constant is -ideal by construction. -""" -from __future__ import annotations - -import math - -from anc.rules import CosineSimilarity, Fixed, Measured -from simroot.pure import PURE_SIMS -from stage20_prediction import MU_GRID, PROFILES, SIGNALS, score - -# the nine functions of the paper, applied to the dissimilarity d = 1 - |cos|, -# plus the identity S(d) = 1 - d which is the raw cosine -FAMILY = [("1 - d (raw cosine)", None, 1.0)] -_SLOPE = {"ratio": 1.0, "exp_abs": 1.0, "rational": 1.0, "ratio_soft": 2.0, - "sqrt": math.inf, "cbrt": math.inf, "rat_sq": 0.0, "exp_sq": 0.0, - "cosine": 0.0} -for _name in ("ratio", "exp_abs", "rational", "ratio_soft", "sqrt", "cbrt", - "rat_sq", "exp_sq", "cosine"): - FAMILY.append((f"S = {_name}", PURE_SIMS[_name], _SLOPE[_name])) - - -def make(S): - r = CosineSimilarity() - r.S = S - return r - - -def main(): - print(__doc__.split("THE BAR")[0]) - - # ---------------- part 1: does the law order this family too? ---------- - print("=" * 88) - print("1. THE FAMILY, ranked - the law of the first domain predicts the " - "order") - print("=" * 88) - print(f"{'similarity':22s}{'slope of (1-S) at 0':>21s}" - f"{'mean NR over 12 cells':>24s}") - print("-" * 88) - scores = [] - for label, S, slope in FAMILY: - vals = [] - for _sig, gen in SIGNALS: - for _pn, prof, _k in PROFILES: - v = score(make(S), gen, prof) - if math.isfinite(v): - vals.append(v) - m = sum(vals) / len(vals) - scores.append((label, slope, m)) - print(f"{label:22s}{slope:>21.3g}{m:>24.2f}", flush=True) - print("-" * 88) - ranked = sorted(scores, key=lambda t: t[2], reverse=True) - print("\n ranked best first:") - for label, slope, m in ranked: - print(f" {label:22s} slope {slope:>8.3g} {m:6.2f} dB") - lin = [m for _l, s, m in scores if s == 1.0] - quad = [m for _l, s, m in scores if s == 0.0] - print(f"\n slope 1 mean {sum(lin) / len(lin):.2f} dB, " - f"slope 0 mean {sum(quad) / len(quad):.2f} dB => the law " - f"{'HOLDS' if min(lin) > max(quad) else 'does NOT hold'} here") - - # ---------------- part 2: the best member against the oracle ----------- - best_label, _s, _m = ranked[0] - best_S = dict((l, S) for l, S, _sl in FAMILY)[best_label] - print("\n" + "=" * 88) - print(f"2. {best_label} against the per-scenario oracle constant") - print("=" * 88) - print(f"{'signal':11s}{'noise profile':17s}{'kind':12s}" - f"{'oracle mu':>11s}{'best mu':>9s}{'variance':>10s}" - f"{'similarity':>12s}{'verdict':>10s}") - print("-" * 88) - wins = {"control": 0, "prediction": 0} - cells = {"control": 0, "prediction": 0} - margins = [] - for signame, gen in SIGNALS: - for pname, prof, kind in PROFILES: - best, arg = -1e30, None - for mu in MU_GRID: - v = score(Fixed(mu), gen, prof) - if math.isfinite(v) and v > best: - best, arg = v, mu - var = score(Measured(), gen, prof) - sim = score(make(best_S), gen, prof) - beats = sim > best - cells[kind] += 1 - wins[kind] += beats - margins.append(sim - best) - print(f"{signame:11s}{pname:17s}{kind:12s}{best:>11.2f}" - f"{arg:>9.2f}{var:>10.2f}{sim:>12.2f}" - f"{'SIMILARITY' if beats else 'oracle':>10s}", flush=True) - tot_w = wins["control"] + wins["prediction"] - tot_c = cells["control"] + cells["prediction"] - print("-" * 88) - print(f"\n beats the per-scenario oracle on {tot_w} of {tot_c} cells" - f" (stationary {wins['control']}/{cells['control']}, " - f"varying {wins['prediction']}/{cells['prediction']})") - print(f" mean margin {sum(margins) / len(margins):+.2f} dB") - print(f"\n the bar was a majority of the twelve. " - f"{'MET' if tot_w > tot_c / 2 else 'NOT MET'}") - - -if __name__ == "__main__": - main() diff --git a/stage22_untuned.py b/stage22_untuned.py deleted file mode 100644 index 1df40a4..0000000 --- a/stage22_untuned.py +++ /dev/null @@ -1,99 +0,0 @@ -"""Stage 22 - the comparison a practitioner actually faces. - -Every comparison so far gave the constant an oracle: the best mu of a grid, -chosen with knowledge of the scenario and the noise profile. That is the -right bar for "can measurement beat perfect tuning?", and the answer is no, -0 of 12. - -But reporting only that bar is its own distortion, the mirror image of -reporting only a bar we pass. A deployable controller does not know the -scenario. The question a practitioner faces is: I must pick mu once, for a -device that will meet conditions I have not seen - how does a fixed choice -compare with measuring? - -So this file asks the same twelve cells a different, and equally -pre-registered, question: - - a single mu, fixed across ALL cells (as a shipped device would have), - against the similarity rule, which is fixed across all cells too because - it has nothing to fix. - -The constant is still given an advantage: it is the best single value over -the whole grid, chosen with hindsight over all twelve cells - a device maker -who guessed perfectly on the first try. What it is not given is a different -value per scenario. - -Both bars are reported together. Neither alone is the answer. -""" -from __future__ import annotations - -import math - -from anc.rules import CosineSimilarity, Fixed, KwongJohnston, Measured -from simroot.pure import PURE_SIMS -from stage20_prediction import MU_GRID, PROFILES, SIGNALS, score - -CELLS = [(sig, gen, pn, prof, kind) - for sig, gen in SIGNALS for pn, prof, kind in PROFILES] - - -def member(name): - r = CosineSimilarity() - r.S = None if name is None else PURE_SIMS[name] - return r - - -def per_cell(rule): - return [score(rule, gen, prof) for _s, gen, _pn, prof, _k in CELLS] - - -def main(): - print(__doc__.split("Both bars")[0]) - print("=" * 90) - - sim = per_cell(member("cbrt")) - var = per_cell(Measured()) - kj = per_cell(KwongJohnston()) - fixed = {mu: per_cell(Fixed(mu)) for mu in MU_GRID} - - # bar 1: the per-cell oracle - a different constant for every scenario - oracle_cell = [max(fixed[mu][i] for mu in MU_GRID) for i in range(len(CELLS))] - # bar 2: the best SINGLE constant across all cells, chosen with hindsight - best_single_mu = max(MU_GRID, key=lambda mu: sum(fixed[mu])) - single = fixed[best_single_mu] - - print(f"{'signal':11s}{'profile':17s}{'per-cell oracle':>17s}" - f"{'best single mu':>16s}{'similarity':>12s}{'KJ':>8s}") - print("-" * 90) - for i, (s, _g, pn, _p, _k) in enumerate(CELLS): - print(f"{s:11s}{pn:17s}{oracle_cell[i]:>17.2f}{single[i]:>16.2f}" - f"{sim[i]:>12.2f}{kj[i]:>8.2f}", flush=True) - print("-" * 90) - n = len(CELLS) - print(f"{'mean':28s}{sum(oracle_cell) / n:>17.2f}{sum(single) / n:>16.2f}" - f"{sum(sim) / n:>12.2f}{sum(kj) / n:>8.2f}") - print(f"\n the best single constant over all twelve cells is " - f"mu = {best_single_mu:g}") - - w_oracle = sum(1 for i in range(n) if sim[i] > oracle_cell[i]) - w_single = sum(1 for i in range(n) if sim[i] > single[i]) - w_kj = sum(1 for i in range(n) if sim[i] > kj[i]) - w_var = sum(1 for i in range(n) if sim[i] > var[i]) - print(f"\n similarity beats the per-cell oracle on {w_oracle:>2d}/{n}" - f" (mean margin {sum(sim[i] - oracle_cell[i] for i in range(n)) / n:+.2f} dB)") - print(f" similarity beats the best single constant on {w_single:>2d}/{n}" - f" (mean margin {sum(sim[i] - single[i] for i in range(n)) / n:+.2f} dB)") - print(f" similarity beats Kwong-Johnston on {w_kj:>2d}/{n}" - f" (mean margin {sum(sim[i] - kj[i] for i in range(n)) / n:+.2f} dB)") - print(f" similarity beats the variance rule on {w_var:>2d}/{n}" - f" (mean margin {sum(sim[i] - var[i] for i in range(n)) / n:+.2f} dB)") - - print("\n Every arbitrary constant, cell by cell:") - print(f" {'mu':>6s}{'mean dB':>10s}{'cells where similarity wins':>30s}") - for mu in MU_GRID: - w = sum(1 for i in range(n) if sim[i] > fixed[mu][i]) - print(f" {mu:>6g}{sum(fixed[mu]) / n:>10.2f}{w:>25d}/{n}") - - -if __name__ == "__main__": - main() diff --git a/stage23_relaxation.py b/stage23_relaxation.py deleted file mode 100644 index c168f0e..0000000 --- a/stage23_relaxation.py +++ /dev/null @@ -1,145 +0,0 @@ -"""Stage 23 - the relaxation parameter, which is where this principle belongs. - -The step rule was stuck on equations that bisection also solves. The wider -statement is: wherever an iterative method has a step length, and there is a -current result and a required one, the step should follow how SIMILAR the two -are. The oldest hand-tuned step length in numerical computing is the -relaxation factor omega of SOR. - - convergence needs 0 < omega < 2 - the optimum known in closed form only for model problems - everywhere else guessed, swept, or left at folklore values - -The similarity form. SOR over-relaxes because the Gauss-Seidel correction -under-shoots consistently in the same direction. How consistent that -direction is, is the cosine between successive corrections, and - - omega = 1 + max(0, cos(d_k, d_{k-1})) lands exactly in [1, 2] - -so the range of the similarity and the admissible interval of the theory -coincide. No constant is chosen. - -The instrument is checked first: on the model Poisson problem the SOR spectral -radius is exactly omega* - 1, so the measured contraction of the oracle must -reproduce it. If it does not, the sweep is wrong and nothing below counts. - -Then three problem classes, in increasing order of how unavailable the optimum -is: - - model Poisson omega* is a theorem - we should LOSE - anisotropic omega* has no closed form - open - variable coefficient / convection - no theory at all - open -""" -from __future__ import annotations - -import math - -import numpy as np - -from relax.sor import (Fixed, OptimalPoisson, SimilarityOmega, poisson_rhs, - residual, solve, sweep_gs) - -TOL = 1e-8 - - -def theory_check(n=16): - """The SOR spectral radius for the model problem is omega* - 1.""" - w = 2.0 / (1.0 + math.sin(math.pi / (n + 1))) - f = poisson_rhs(n) - h2 = 1.0 / (n + 1) ** 2 - u = np.zeros_like(f) - rs = [] - for _ in range(60): - sweep_gs(u, f, h2, w) - rs.append(float(np.linalg.norm(residual(u, f, h2)))) - ratios = [rs[i + 1] / rs[i] for i in range(30, len(rs) - 1) if rs[i] > 0] - measured = sum(ratios) / len(ratios) if ratios else float("nan") - return w, w - 1.0, measured - - -def anisotropic(n, eps): - """-u_xx - eps*u_yy: the classic problem where a single omega is wrong.""" - cw = np.full((n, n), 1.0) - ce = np.full((n, n), 1.0) - cs = np.full((n, n), eps) - cn = np.full((n, n), eps) - cc = np.full((n, n), 2.0 + 2.0 * eps) - return (cw, ce, cs, cn, cc) - - -def variable(n, rng): - """A jumping diffusion coefficient - no closed-form optimum exists.""" - a = np.exp(rng.uniform(-2.0, 2.0, (n, n))) - cw = np.roll(a, 1, 0) * 0.5 + a * 0.5 - ce = np.roll(a, -1, 0) * 0.5 + a * 0.5 - cs = np.roll(a, 1, 1) * 0.5 + a * 0.5 - cn = np.roll(a, -1, 1) * 0.5 + a * 0.5 - cc = cw + ce + cs + cn - return (cw, ce, cs, cn, cc) - - -def convection(n, pe): - """Convection-diffusion, upwind: non-symmetric, and the SOR theory for - symmetric problems does not apply at all.""" - cw = np.full((n, n), 1.0 + pe) - ce = np.full((n, n), 1.0) - cs = np.full((n, n), 1.0) - cn = np.full((n, n), 1.0) - cc = np.full((n, n), 4.0 + pe) - return (cw, ce, cs, cn, cc) - - -def main(): - n = 24 - w_opt, rho_theory, rho_measured = theory_check(n) - print("instrument check on the model Poisson problem") - print(f" omega* = {w_opt:.6f} spectral radius: theory " - f"{rho_theory:.6f}, measured {rho_measured:.6f}") - if abs(rho_theory - rho_measured) > 0.02: - print(" the sweep does NOT reproduce the theory - stopping") - return - print(" the sweep reproduces the theory\n") - - rng = np.random.default_rng(7) - f = poisson_rhs(n, rng) - grid = [ - ("model Poisson", None, "omega* is a theorem"), - ("anisotropic eps=0.1", anisotropic(n, 0.1), "no closed form"), - ("anisotropic eps=0.01", anisotropic(n, 0.01), "no closed form"), - ("variable coefficient", variable(n, np.random.default_rng(3)), - "no theory"), - ("convection Pe=2", convection(n, 2.0), "non-symmetric"), - ("convection Pe=10", convection(n, 10.0), "non-symmetric"), - ] - sweep_grid = (1.0, 1.2, 1.4, 1.6, 1.8, 1.9, w_opt) - - print("=" * 96) - print("sweeps to ||r|| <= 1e-8 ||r0|| (lower is better; 20000 = did " - "not converge)") - print("=" * 96) - print(f"{'problem':24s}{'why':16s}{'Gauss-Seidel':>14s}" - f"{'best fixed w':>14s}{'(that w)':>10s}{'similarity':>12s}" - f"{'mean w':>8s}") - print("-" * 96) - wins = total = 0 - for label, coeff, why in grid: - gs, _ = solve(f, Fixed(1.0), tol=TOL, coeff=coeff) - best, bw = 10 ** 9, None - for w in sweep_grid: - k, _ = solve(f, Fixed(w), tol=TOL, coeff=coeff) - if k < best: - best, bw = k, w - sim, mw = solve(f, SimilarityOmega(), tol=TOL, coeff=coeff) - total += 1 - wins += sim < best - print(f"{label:24s}{why:16s}{gs:>14d}{best:>14d}{bw:>10.3f}" - f"{sim:>12d}{mw:>8.3f}", flush=True) - print("-" * 96) - print(f"\n similarity beats the best swept omega on {wins} of {total} " - f"problems") - print(" the swept omega is an oracle: it costs seven full solves to " - "find,\n and the similarity rule costs none.") - - -if __name__ == "__main__": - main() diff --git a/tests/test_anc.py b/tests/test_anc.py deleted file mode 100644 index 8c9ef83..0000000 --- a/tests/test_anc.py +++ /dev/null @@ -1,91 +0,0 @@ -"""Active noise control: the step size of an NLMS canceller. - -The third domain, and the third instance of the same pattern. What is pinned -here is the pattern and the mistake, not a victory. -""" -from __future__ import annotations - -import math -import os -import sys - -import numpy as np -import pytest - -sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) - -from anc.bench import SCENARIOS, make_primary_path, make_secondary_path, run_nlms -from anc.rules import Fixed, InSampleMeasured, KwongJohnston, Measured - -N, SEEDS = 8000, 4 - - -SNRS = (30.0, 15.0) - - -def sweep(rule, snrs=SNRS): - """Both noise levels, as the experiment uses. - - Measuring at 30 dB alone is the wrong test and was: with almost no - measurement noise there is nothing for a safety factor to hedge against, - the exact step mu = 1 is near optimal, and any hedging costs a little. - The whole point of the rule appears at 15 dB. A claim about the grid has - to be tested on the grid.""" - vals = [] - for snr in snrs: - for s in range(SEEDS): - for _sname, gen in SCENARIOS: - rng = np.random.default_rng(1000 + s) - noise = 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), min(vals), len(vals) - - -def test_the_in_sample_observable_is_uninformative(): - """The mistake, pinned. NLMS with mu = 1 drives the a posteriori error of - the current sample to exactly zero, so a rule that measures it reads 'the - model was perfect' and proposes mu = 1 for ever. It is therefore the - fixed mu = 1 rule under another name, and the numbers say so.""" - a, _lo, _n = sweep(InSampleMeasured()) - b, _lo2, _n2 = sweep(Fixed(1.0)) - assert abs(a - b) < 0.1, (a, b) - - -def test_measuring_out_of_sample_is_better_than_measuring_in_sample(): - """The fix: what mu hedges against is the part of the error the filter - cannot cancel, and that is observable as the floor the error keeps - hitting - not as the residual of the sample just fitted.""" - good, _lo, _n = sweep(Measured()) - bad, _lo2, _n2 = sweep(InSampleMeasured()) - assert good > bad, (good, bad) - - -def test_the_measured_rule_beats_the_four_constant_scheme(): - """Kwong-Johnston is the classical variable-step-size rule and carries - four constants. This is the comparison the transfer is about.""" - ours, _lo, _n = sweep(Measured()) - kj, _lo2, _n2 = sweep(KwongJohnston()) - assert ours > kj + 3.0, (ours, kj) - - -def test_and_loses_to_one_well_chosen_constant(): - """The other half, and the same in all three domains: measuring the - constant beats the elaborate multi-constant schemes and does not beat a - single constant chosen correctly for the benchmark.""" - ours, _lo, _n = sweep(Measured()) - best_fixed, _lo2, _n2 = sweep(Fixed(0.5)) - assert best_fixed > ours, (best_fixed, ours) - assert best_fixed - ours < 3.0, (best_fixed, ours) - - -def test_no_rule_diverges(): - """A step size that is too large makes the filter blow up; none of these - may.""" - for rule in (Fixed(0.1), Fixed(0.5), Fixed(1.0), KwongJohnston(), - Measured()): - _m, _lo, n = sweep(rule) - assert n == SEEDS * len(SCENARIOS) * len(SNRS), (rule.name, n) diff --git a/tests/test_bulletproof.py b/tests/test_bulletproof.py index 19b0c16..cc5ba82 100644 --- a/tests/test_bulletproof.py +++ b/tests/test_bulletproof.py @@ -123,39 +123,3 @@ def test_the_scope_of_the_claim_is_repairs_not_the_field(): 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) diff --git a/tests/test_ode_controllers.py b/tests/test_ode_controllers.py deleted file mode 100644 index f59d174..0000000 --- a/tests/test_ode_controllers.py +++ /dev/null @@ -1,102 +0,0 @@ -"""The second domain: adaptive step-size control for ODEs. - -What is being tested is the transfer, not the integrator: everything except -the controller is shared, so a difference can only come from the controller. -""" -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 odestep.controllers import (ALL, Classical, Elementary, GustafssonPI, - MeasuredSafety) -from odestep.driver import integrate, validate_against_scipy -from stage15_ode import PROBLEMS, TOLS, reference, scaled_err - - -def test_the_driver_reproduces_the_library(): - """The instrument before the experiment. - - The first version of this work carried a hand-typed Dormand-Prince - tableau and validated it by measuring convergence orders. That was both - an unnecessary risk and, as it turned out, a material one: with the - hand-written tableau the measured-safety controller appeared to beat the - five-constant PI controller by 8%, and with SciPy's it is 2.6% behind. - The integrator is now SciPy's throughout and the check is direct.""" - for name, mine, ref, endpoint in validate_against_scipy(): - assert abs(mine - ref) / ref < 0.05, (name, mine, ref) - assert endpoint < 1e-6, (name, endpoint) - - -def test_every_controller_solves_an_easy_problem(): - f = lambda t, y: [-y[0]] - for c in ALL: - r = integrate(f, 0.0, 10.0, [1.0], c, atol=1e-8, rtol=1e-8) - assert r["reached_end"], c.name - assert abs(r["y"][0] - math.exp(-10.0)) < 1e-6, c.name - - -def test_the_constant_free_baseline_is_the_one_to_beat(): - """The elementary controller is what the error model alone dictates and - what every practical code refuses to use. If it were not much worse, - there would be nothing for a safety factor - chosen or measured - to buy.""" - f = lambda t, y: [y[1], 5.0 * (1 - y[0] * y[0]) * y[1] - y[0]] - elem = integrate(f, 0.0, 25.0, [2.0, 0.0], Elementary(), atol=1e-8, - rtol=1e-8) - meas = integrate(f, 0.0, 25.0, [2.0, 0.0], MeasuredSafety(), atol=1e-8, - rtol=1e-8) - # margins measured with SciPy's integrator. An earlier version demanded - # a factor of two on evaluations and five on rejections; those came from - # a hand-written tableau and are not what the library produces. - assert elem["nfev"] > 1.2 * meas["nfev"], (elem["nfev"], meas["nfev"]) - assert elem["rejected"] > 1.5 * meas["rejected"], (elem["rejected"], - meas["rejected"]) - - -def test_where_the_measured_safety_factor_actually_stands(): - """The honest position, over the whole grid, with SciPy's integrator. - - It is NOT ahead of the PI controller - an earlier version of this test - asserted that, on numbers produced by a hand-written tableau. What it is: - level with a five-constant controller to within a few percent while - having none, and clearly ahead of the constant-free elementary rule that - every practical code refuses to use. It is behind the three-constant - classical controller and the test says so.""" - tot = {c.name: 0 for c in ALL} - fails = {c.name: 0 for c in ALL} - refs = {name: (exact if exact is not None else reference(f, t0, t1, y0)) - for name, f, t0, t1, y0, exact in PROBLEMS} - for tol in TOLS: - for name, f, t0, t1, y0, _exact in PROBLEMS: - for c in ALL: - r = integrate(f, t0, t1, y0, c, atol=tol, rtol=tol) - err = scaled_err(r["y"], refs[name], tol, tol) - bad = (not r["reached_end"]) or (math.isfinite(err) - and err > 1000.0) - tot[c.name] += r["nfev"] - fails[c.name] += bad - # nobody wins by failing more - assert len(set(fails.values())) == 1, fails - m = tot[MeasuredSafety().name] - pi = tot[GustafssonPI().name] - # level with the five-constant PI controller, within five percent - assert abs(m - pi) / pi < 0.05, (m, pi) - # clearly ahead of the constant-free rule the error model alone dictates - assert m < 0.9 * tot[Elementary().name], (m, tot[Elementary().name]) - # and honestly behind the three-constant classical one - assert m > tot[Classical().name], (m, tot[Classical().name]) - - -def test_the_controller_holds_no_unused_state(): - """It was called a credit controller while banking a margin that took no - part in any decision - the ledger was vestigial and the name unearned. - Nothing may accumulate here that the step factor does not read.""" - c = MeasuredSafety() - c.reset() - fields = set(vars(c)) - assert fields == {"predicted", "delta"}, fields diff --git a/tests/test_registered_prediction.py b/tests/test_registered_prediction.py deleted file mode 100644 index 9372d21..0000000 --- a/tests/test_registered_prediction.py +++ /dev/null @@ -1,77 +0,0 @@ -"""A prediction registered before the run, and falsified by it. - -The hypothesis was that a method which wins in some places and loses in others -must have a regime where it wins outright, and that the regime is where the -quantity the constant approximates CHANGES during the run - because one -constant cannot be right at two different noise levels while a measured value -tracks. - -It was stated before measuring, with the falsification condition attached: the -measured rule had to lose every constant-noise cell and win every -varying-noise cell. It lost all twelve. - -The first version of the test could not have confirmed it either way - noise -reduction was scored on the second half of the run, by which time every ramp -and step had already arrived, so the metric never saw a transition. That flaw -was found and fixed before the conclusion was drawn, and the conclusion did -not change. - -This file exists so the line is not reopened without new evidence. -""" -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 anc.rules import Fixed, Measured -from stage20_prediction import (FLAT_5, FLAT_30, STEP_25, SWING_25, sig_cabin, - sig_tonal, score) - -SEEDS = 3 - - -def test_the_metric_window_sees_the_transition(): - """The flaw that would have produced a false falsification. A varying - profile must not score the same as a constant one at its final level.""" - varying = score(Fixed(0.1), sig_tonal, SWING_25, seeds=SEEDS) - final_level = score(Fixed(0.1), sig_tonal, FLAT_5, seeds=SEEDS) - assert varying > final_level + 2.0, (varying, final_level) - - -@pytest.mark.parametrize("signal", [sig_tonal, sig_cabin], - ids=["tonal", "cabin"]) -@pytest.mark.parametrize("profile,kind", [(FLAT_30, "control"), - (FLAT_5, "control"), - (SWING_25, "prediction"), - (STEP_25, "prediction")], - ids=["flat30", "flat5", "ramp", "step"]) -def test_an_oracle_constant_beats_the_measured_rule_everywhere(signal, profile, - kind): - """The registered prediction required the measured rule to WIN the two - 'prediction' cells. It does not, on either signal. A single constant, - chosen with knowledge of the scenario, is ahead even when the noise level - swings by 25 dB inside the run.""" - ms = score(Measured(), signal, profile, seeds=SEEDS) - best = max(score(Fixed(mu), signal, profile, seeds=SEEDS) - for mu in (0.05, 0.1, 0.2, 0.5, 1.0)) - assert best > ms, (kind, best, ms) - - -def test_the_conclusion_this_licenses(): - """Winning somewhere and losing elsewhere does not imply a regime of - outright victory exists; it is equally consistent with a middling - generalist. On this evidence there is no such regime, and the claim the - work may make is the scoped one in test_modern_repairs.py.""" - losses = 0 - for signal in (sig_tonal, sig_cabin): - for profile in (FLAT_30, FLAT_5, SWING_25, STEP_25): - ms = score(Measured(), signal, profile, seeds=SEEDS) - best = max(score(Fixed(mu), signal, profile, seeds=SEEDS) - for mu in (0.05, 0.2, 1.0)) - losses += ms < best - assert losses == 8, losses diff --git a/tests/test_relaxation.py b/tests/test_relaxation.py deleted file mode 100644 index 7ad9f86..0000000 --- a/tests/test_relaxation.py +++ /dev/null @@ -1,95 +0,0 @@ -"""The relaxation parameter: the fourth domain, and the fourth time the same -pattern appears. - -omega in SOR is the oldest hand-tuned step length in numerical computing. The -similarity form - omega = 1 + cos between successive corrections - lands -exactly in the admissible interval [1, 2] with nothing chosen, and it beats -Gauss-Seidel five-fold. It loses to a swept omega two-fold. - -A second map, derived from the theory rather than from the geometry, was tried -and is worse than the first. Both are pinned: the working one so it is not -lost, the refuted one so it is not retried. -""" -from __future__ import annotations - -import math -import os -import sys - -import numpy as np -import pytest - -sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) - -from relax.sor import (ContractionOmega, Fixed, SimilarityOmega, poisson_rhs, - residual, solve, sweep_gs) -from stage23_relaxation import anisotropic, convection, theory_check - -N = 14 -F = poisson_rhs(N) -W_OPT = 2.0 / (1.0 + math.sin(math.pi / (N + 1))) - - -def test_the_sweep_reproduces_the_sor_theory(): - """The instrument before the experiment: for the model problem the SOR - spectral radius is exactly omega* - 1.""" - _w, theory, measured = theory_check(16) - assert abs(theory - measured) < 0.02, (theory, measured) - - -def test_the_similarity_omega_lands_in_the_admissible_interval(): - """Range of the similarity and convergence condition of the method - coincide, which is why nothing had to be clipped by hand.""" - r = SimilarityOmega() - r.reset() - d1 = np.array([1.0, 2.0, 3.0]) - assert r.next_omega(None, None) == 1.0 - assert r.next_omega(d1, d1) == pytest.approx(2.0) # identical - assert r.next_omega(d1, -d1) == pytest.approx(1.0) # opposed - assert r.next_omega(np.array([1.0, 0.0, 0.0]), - np.array([0.0, 1.0, 0.0])) == pytest.approx(1.0) - - -@pytest.mark.parametrize("label,coeff", [ - ("poisson", None), ("aniso", anisotropic(N, 0.1))]) -def test_similarity_beats_gauss_seidel_by_a_wide_margin(label, coeff): - gs, _ = solve(F, Fixed(1.0), coeff=coeff) - sim, _ = solve(F, SimilarityOmega(), coeff=coeff) - assert sim * 3 < gs, (label, sim, gs) - - -@pytest.mark.parametrize("label,coeff", [ - ("poisson", None), ("aniso", anisotropic(N, 0.1))]) -def test_and_loses_to_a_swept_omega(label, coeff): - """The other half, stated so the claim is never inflated.""" - sim, _ = solve(F, SimilarityOmega(), coeff=coeff) - best = min(solve(F, Fixed(w), coeff=coeff)[0] - for w in (1.2, 1.4, 1.6, 1.8, W_OPT)) - assert best < sim, (label, best, sim) - - -def test_the_similarity_rule_over_relaxes(): - """The diagnosis, not just the outcome: the cosine between successive SOR - corrections sits near 1, so 1 + cos saturates near 2 and the rule picks a - larger omega than the optimum.""" - _sim, mean_w = solve(F, SimilarityOmega()) - assert mean_w > W_OPT, (mean_w, W_OPT) - - -def test_the_contraction_map_is_refuted_and_worse(): - """The second attempt, pinned so it is not tried again. Undamped - feedback pins omega at the limit from the first sweep.""" - sim, _ = solve(F, SimilarityOmega()) - con, mean_w = solve(F, ContractionOmega()) - assert con > sim, (con, sim) - assert mean_w > 1.9, mean_w - - -def test_convection_is_where_it_actually_hurts(): - """On a strongly non-symmetric problem Gauss-Seidel is already fast and - the rule makes it slower - the one case where the similarity form is not - merely second but harmful.""" - coeff = convection(N, 10.0) - gs, _ = solve(F, Fixed(1.0), coeff=coeff) - sim, _ = solve(F, SimilarityOmega(), coeff=coeff) - assert sim > gs, (sim, gs) diff --git a/tests/test_similarity_family_anc.py b/tests/test_similarity_family_anc.py deleted file mode 100644 index f230f44..0000000 --- a/tests/test_similarity_family_anc.py +++ /dev/null @@ -1,105 +0,0 @@ -"""The similarity family for the NLMS step size, and what it settled. - -Two corrections came from outside and both were right: - - 1. mu should be a SIMILARITY, not a variance estimate. The filter can only - remove what correlates with its input, so the quantity is the cosine - between the error and the filtered reference. It beats the variance - formulation on eleven of twelve cells. - - 2. Testing one similarity function is as incomplete here as it would have - been in the root-finding work. The whole family must be run, and when it - is, the best member is not the cosine. - -And the law of the first domain - that the slope of (1 - S) at the origin -decides, not the shape of S - orders this family too, monotonically and -without exception. That is a prediction made in one domain and confirmed in -another. - -What it did not settle: even the best member loses to a constant tuned per -scenario, on all twelve cells. -""" -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 anc.rules import CosineSimilarity, Fixed, Measured -from simroot.pure import PURE_SIMS -from stage20_prediction import FLAT_15, FLAT_30, STEP_25, sig_cabin, sig_tonal, \ - score - -SEEDS = 3 -CELLS = [(sig_tonal, FLAT_30), (sig_tonal, STEP_25), - (sig_cabin, FLAT_15), (sig_cabin, STEP_25)] - - -def make(name): - r = CosineSimilarity() - r.S = None if name is None else PURE_SIMS[name] - return r - - -def mean(rule): - vals = [score(rule, gen, prof, seeds=SEEDS) for gen, prof in CELLS] - return sum(vals) / len(vals) - - -def test_the_similarity_formulation_beats_the_variance_one(): - """Correction 1. mu as 'how much of the error resembles what the filter - can remove' against mu as 'one minus the noise-power fraction'. - - Stated carefully, because the first version of this test was not: what - beats the variance rule on eleven of twelve cells is the best member of - the FAMILY, not the raw cosine. The raw cosine can lose to it - on this - subset it does, 14.04 against 14.32 - which is precisely why correction 2 - (run the whole family) was needed before the claim could be made at all. - """ - best_member = max(mean(make(n)) for n in (None, "ratio", "exp_abs", - "rational", "ratio_soft", - "sqrt", "cbrt")) - assert best_member > mean(Measured()), best_member - - -def test_the_raw_cosine_is_not_the_best_member(): - """Correction 2. A single member is not the family: the identity - S(d) = 1 - d is beaten by members with a steeper approach to 1.""" - raw = mean(make(None)) - best = max(mean(make(n)) for n in ("sqrt", "cbrt", "ratio_soft")) - assert best > raw, (best, raw) - - -def test_two_code_paths_for_the_same_function_agree(): - """S = None takes mu = |cos| directly; S = ratio computes - max(0, 1 - d) with d = 1 - |cos|. They are the same function on [0,1] - and reach it by different code, so they must score identically.""" - a = mean(make(None)) - b = mean(make("ratio")) - assert abs(a - b) < 1e-9, (a, b) - - -def test_the_law_of_the_first_domain_orders_this_family_too(): - """The prediction carried across domains: slope of (1 - S) at 0 decides. - Every slope-1 member must beat every slope-0 member.""" - slope1 = [mean(make(n)) for n in ("ratio", "exp_abs", "rational")] - slope0 = [mean(make(n)) for n in ("rat_sq", "exp_sq", "cosine")] - assert min(slope1) > max(slope0), (slope1, slope0) - - -def test_and_it_still_loses_to_a_tuned_constant(): - """The half that did not work, pinned so the claim is never overstated. - The competitor knows the scenario; no deployable controller does - and it - is still ahead.""" - losses = 0 - for gen, prof in CELLS: - best_member = max(score(make(n), gen, prof, seeds=SEEDS) - for n in (None, "cbrt", "sqrt", "ratio")) - oracle = max(score(Fixed(mu), gen, prof, seeds=SEEDS) - for mu in (0.05, 0.1, 0.2, 0.5, 1.0)) - losses += best_member < oracle - assert losses == len(CELLS), losses