Compare commits
1
Commits
886b3f4244
..
main
| Author | SHA1 | Date | |
|---|---|---|---|
  vvsandClaude Opus 5
|
8ad09b47a6 |
@@ -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