The three neighbouring applications out of the published branch: the NLMS step size (anc), the ODE step controllers (odestep) and the relaxation factor (relax), with the six stage scripts that drive them and the five tests that cover them. They are not part of the root-finding paper and no document refers to them. stage18_bulletproof.py keeps claim A and loses claim B; tests/test_bulletproof.py loses the one test that reached into anc. Server: verify 548 OK / 0 BAD, tests x2 1205 passed (was 1240 with the three).
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -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
|
||||
|
||||
-190
@@ -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)
|
||||
-204
@@ -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>| / (||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 <e, x>
|
||||
self.ee = None # running <e, e>
|
||||
self.xx = None # running <x, x>
|
||||
|
||||
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()]
|
||||
@@ -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()]
|
||||
@@ -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
|
||||
-238
@@ -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))
|
||||
-189
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
|
||||
@@ -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()
|
||||
@@ -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>| / (||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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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)
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user