vvsandClaude Opus 5
2357379c2d
Not empty-handed: our own parameter-free reading 'defect' - refuse when |phi - xi| > min(xi, 1 - xi) - sets a new canonical best, APS 2523 at 153/154 (previous 2527, Chandrupatla 2785), guarantee intact (0 of 135 family cells), refusing on 7.0% of steps. It is a second point on the same trade, not a replacement: 5733 over both sets against 5511, because CHA gives back 3210. Also: _shape_ratios crashed on a missing third point, which evaluate would have turned into a silent sweep of zeros; the slope-gate pin was anchored on line spelling; and the server copy had drifted (stale relax/sor.py, test_relaxation.py absent - 746 tests instead of 755). All synced, 92/92 files identical by content, 755 pass on the server. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
918 lines
48 KiB
Python
918 lines
48 KiB
Python
"""The pure similarity root finder: one similarity function, one step rule.
|
|
|
|
direction : the side of the bracket that contains the sign change
|
|
length : (1 - S(z)) * (distance to the far end), z = f(x) / scale
|
|
|
|
Everything about the method is in S and in where `scale` comes from.
|
|
|
|
scale CHOSEN (a fixed reference value) -> linear convergence at a rate
|
|
that is a property of the problem, not of the method, and
|
|
therefore fails on some problems (stage 0b, 5-8 of 8)
|
|
scale DERIVED (the local slope of f, read off the two most recent
|
|
evaluations) -> the step becomes |f|/|f'| and
|
|
the convergence is superlinear (stage 0d/0e)
|
|
|
|
The derivation, from the contraction law measured in stage 0c:
|
|
|
|
e_{n+1} = e_n * |1 - c*|f'(r)|*D/scale| , c = slope of (1 - S) at 0
|
|
rate = 0 <=> c = scale / (|f'(r)|*D)
|
|
|
|
so the one quantity the pure rule needs is the LOCAL SCALE OF THE VALUE AXIS.
|
|
Choosing it is what makes the method brittle; reading it off the data is what
|
|
makes it parameter free - the same move the SIM control paper makes when it
|
|
normalises the input space from the data instead of choosing membership
|
|
widths by hand.
|
|
|
|
Safeguards
|
|
----------
|
|
An interpolation-only rule has no worst-case bound. Two parameter-free
|
|
safeguards are provided and are measured against each other:
|
|
|
|
"bracket" bisect whenever a similarity step failed to halve the BRACKET.
|
|
Sound, but expensive: a secant step routinely collapses the
|
|
error while barely moving the proven bracket, so the guard
|
|
fires on steps that were in fact excellent (measured: 150 vs
|
|
86 evaluations over the stage-0 problem set).
|
|
|
|
"step" bisect whenever the proposed step is not at most half the step
|
|
taken two evaluations ago. Step lengths are then forced to
|
|
decrease geometrically, and a converging secant sequence -
|
|
whose steps shrink far faster than by half - never triggers
|
|
it. This is the criterion Brent's method uses internally,
|
|
stated without his tuning constants.
|
|
|
|
It does NOT on its own bound the worst case, and an earlier
|
|
version of this docstring claimed that it did. A sequence
|
|
that creeps towards a non-root from one side with steps that
|
|
halve every time has small steps by construction, so the guard
|
|
stays silent, the sign never changes, and the proven bracket
|
|
never moves. Cheap, but not sound alone.
|
|
|
|
"both" the step guard, with the bracket guard underneath it as the
|
|
backstop that actually carries the worst-case bound. The step
|
|
guard does the work on ordinary problems, so the bracket guard
|
|
almost never fires and almost never costs anything.
|
|
|
|
"project" do not REPLACE the similarity step at all - PROJECT it onto the
|
|
largest interval still consistent with the bisection bound:
|
|
|
|
allowed_j = 2^(n_max - j) * (xtol/2)
|
|
x = clamp(x_sim, hi - allowed_j, lo + allowed_j)
|
|
|
|
with n_max = 2*ceil(log2(w0 / xtol)) fixed by the initial
|
|
bracket. The next bracket is at most max(x - lo, hi - x), so
|
|
the clamp gives w_{j+1} <= allowed_j directly, and induction
|
|
bounds the run by n_max steps in exact arithmetic - one more
|
|
when the bracket lands exactly on the boundary in floating
|
|
point. The guarantee is therefore
|
|
2*ceil(log2(w0/xtol)) + 3 evaluations: twice bisection plus
|
|
one step.
|
|
|
|
Zero slack (n_max = ceil(log2(w0/xtol)), a worst case of
|
|
exactly bisection) is degenerate and was measured as such: it
|
|
demands a halving at every single step, leaving the projection
|
|
no freedom, and the guard behaves as bisection.
|
|
|
|
The three guards above discard a computed step; this one keeps
|
|
as much of it as the bound allows, which is why it does not
|
|
flatten the order. The construction is ITP's projection
|
|
stage. ITP's three constants (kappa1, kappa2, n0) all live in
|
|
its TRUNCATION stage, which exists to repair regula falsi;
|
|
the similarity step needs no repair, so the truncation - and
|
|
with it every constant - is dropped.
|
|
|
|
"safe" the step guard for cheapness, with the projection carrying the
|
|
bound. This is the configuration the paper proposes: a proven
|
|
worst case AND the full order, which the three earlier guards
|
|
could not deliver together.
|
|
|
|
The only constant in any of them is 1/2, which is the definition of bisection
|
|
rather than a tuned value.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from .baselines import Counter, _Budget
|
|
|
|
# --------------------------------------------------------------------------
|
|
# [0,1] similarity functions of a normalised residual z: S(0) = 1 (identical
|
|
# to the expectation), S -> 0 far from it.
|
|
# --------------------------------------------------------------------------
|
|
PURE_SIMS = {
|
|
# the ratio form S = 1 - |f|/scale. Note that the symmetric ratio
|
|
# min(|x|,|y|)/max(|x|,|y|) is identically 0 when one argument is the
|
|
# target, so this is the only usable reading of "ratio similarity to 0".
|
|
"ratio": lambda z: max(0.0, 1.0 - abs(z)),
|
|
"sqrt": lambda z: max(0.0, 1.0 - math.sqrt(abs(z))),
|
|
"cbrt": lambda z: max(0.0, 1.0 - abs(z) ** (1.0 / 3.0)),
|
|
"rational": lambda z: 1.0 / (1.0 + abs(z)),
|
|
"ratio_soft": lambda z: 1.0 / (1.0 + abs(z)) ** 2,
|
|
"exp_abs": lambda z: math.exp(-abs(z)),
|
|
# similarities that approach 1 quadratically: c = 0, predicted rate 1,
|
|
# kept because the law has to be falsifiable in both directions
|
|
"rat_sq": lambda z: 1.0 / (1.0 + z * z),
|
|
"exp_sq": lambda z: math.exp(-z * z),
|
|
"cosine": lambda z: 1.0 / math.sqrt(1.0 + z * z),
|
|
}
|
|
|
|
GUARDS = ("none", "bracket", "step", "both", "project", "safe", "credit")
|
|
|
|
SPACES = ("x", "log")
|
|
|
|
|
|
POINT_RULES = ("recent", "shape")
|
|
DIRECTIONS = ("distance", "dissimilar")
|
|
|
|
EXPONENTS = ("none", "order")
|
|
|
|
|
|
def _read_z(resid, sc, exponent, m):
|
|
"""The normalised residual the similarity function is applied to.
|
|
|
|
Valeri, 2026-09-21: "the similarity function may simply be the wrong one", and then
|
|
"the similarity function can differ in different PLACES - not by domain, by the place
|
|
of the parameter". This is the step-length place, and here the reading is fixed by
|
|
what f looks like locally rather than chosen.
|
|
|
|
For |f| ~ K.e^m the residual ratio against the far end is dimensionless and needs no
|
|
slope at all: z = |f(x)|/|f(far)| = (e/E)^m, so z^(1/m) = e/E recovers the DISTANCE
|
|
ratio exactly. Under S = 1/(1+z) the step fraction is z/(1+z), which at m = 1 is
|
|
false position (Proposition 4 in the paper: 127/154, 10560 evaluations) - and false
|
|
position's three-century-old stall is precisely a similarity read at order 1 on an
|
|
order-m approach. Reading it at the MEASURED order turns the same formula into
|
|
order-m inverse interpolation, exact in one step for a pure power.
|
|
|
|
m comes from _local_order, i.e. from three evaluations already paid for, and is 1
|
|
whenever the power model does not apply, so nothing here is selected.
|
|
"""
|
|
z = resid / sc
|
|
if exponent == "order" and m is not None and m > 1.0 and z > 0.0:
|
|
try:
|
|
z = z ** (1.0 / m)
|
|
except (OverflowError, ValueError):
|
|
return resid / sc
|
|
if not math.isfinite(z):
|
|
return resid / sc
|
|
return z
|
|
|
|
|
|
def _shape_ratios(x, fx, u, fu, v, fv):
|
|
"""(xi, phi) for the triple, or None when it is degenerate.
|
|
|
|
The two dimensionless ratios the DECISION place is written in (PREREG 2026-09-21,
|
|
slot 2): xi of positions, phi of values, equal for a linear f. _shape_defect is
|
|
|phi - xi| of these same two; both callers must read one definition, or the step
|
|
place and the decision place could silently disagree about what the triple is.
|
|
|
|
A missing point (the first steps of a run have no third evaluation yet) is a
|
|
degenerate triple, not an error: returning None here is what keeps every caller from
|
|
raising inside the loop, where `evaluate` would swallow it and charge the full budget
|
|
- the failure mode that produced a whole sweep of 39800 on 2026-09-20."""
|
|
if x is None or u is None or v is None:
|
|
return None
|
|
if v == u or fv == fu:
|
|
return None
|
|
xi = (x - u) / (v - u)
|
|
ph = (fx - fu) / (fv - fu)
|
|
if not (math.isfinite(xi) and math.isfinite(ph)):
|
|
return None
|
|
return xi, ph
|
|
|
|
|
|
def _decide_refuse(mode, S, z, x, fx, u, fu, v, fv):
|
|
"""The DECISION place: does the local model deserve this evaluation, or is halving
|
|
the better use of it? Returns True to refuse the similarity step.
|
|
|
|
Three readings, all of quantities already in hand and all without a constant:
|
|
|
|
"band" Chandrupatla's 1997 acceptance band, 1 - sqrt(1-xi) < phi < sqrt(xi).
|
|
NOT OURS. It is in the grid as an INSTRUMENT: it answers whether the
|
|
remaining gap lives in this place at all. It must never be presented as
|
|
part of the proposed method.
|
|
"defect" refuse when the configuration is less like a linear one than the iterate
|
|
is off-centre: |phi - xi| > min(xi, 1 - xi). Two computed dimensionless
|
|
quantities compared; no threshold.
|
|
"simpair" refuse when the configuration's similarity to linear is below the
|
|
iterate's similarity to the root: S(|phi - xi|) < S(z). The SAME S is
|
|
used in both places on purpose - pairing a different S per place would be
|
|
a menu, and a menu is the same fault as a hardcoded constant.
|
|
|
|
A degenerate triple carries no decision, so the step is allowed and the safeguard
|
|
keeps it honest; refusing on missing information would be a guess.
|
|
"""
|
|
r = _shape_ratios(x, fx, u, fu, v, fv)
|
|
if r is None:
|
|
return False
|
|
xi, ph = r
|
|
if mode == "band":
|
|
if not (0.0 < xi < 1.0):
|
|
return True
|
|
return not (1.0 - math.sqrt(1.0 - xi) < ph < math.sqrt(xi))
|
|
defect = abs(ph - xi)
|
|
if mode == "defect":
|
|
return defect > min(abs(xi), abs(1.0 - xi))
|
|
if mode == "simpair":
|
|
return S(defect) < S(z)
|
|
return False
|
|
|
|
|
|
def _shape_defect(x, fx, u, fu, v, fv):
|
|
"""|phi - xi| for the triple, or None when it is degenerate.
|
|
|
|
xi = (x - u)/(v - u) is a ratio of POSITIONS, phi = (fx - fu)/(fv - fu) a ratio of
|
|
VALUES; both are invariant under affine rescaling of x and of f. For a linear f they
|
|
are equal, so their gap measures how unlike a linear configuration the triple is —
|
|
the same dimensionless currency Chandrupatla's acceptance band is written in."""
|
|
r = _shape_ratios(x, fx, u, fu, v, fv)
|
|
return None if r is None else abs(r[1] - r[0])
|
|
|
|
|
|
def _pick_triple(x, fx, x_prev, f_prev, x_prev2, f_prev2, A, FA, B, FB, far):
|
|
"""SELECT the triple whose configuration is least unlike a linear one (PREREG
|
|
2026-09-21, mode A). Candidates are declared, not searched: the recency triple, the
|
|
recent point with the far bracket end, and the two proven bracket ends. Returns
|
|
((xu, fu), (xv, fv)) or None when nothing is admissible."""
|
|
far_f = FB if far == B else FA
|
|
cands = (((x_prev, f_prev), (x_prev2, f_prev2)),
|
|
((x_prev, f_prev), (far, far_f)),
|
|
((A, FA), (B, FB)))
|
|
best, best_d = None, None
|
|
for u, v in cands:
|
|
if u[0] is None or v[0] is None:
|
|
continue
|
|
d = _shape_defect(x, fx, u[0], u[1], v[0], v[1])
|
|
if d is None or x == u[0] or x == v[0]:
|
|
continue
|
|
if best_d is None or d < best_d:
|
|
best, best_d = (u, v), d
|
|
return best
|
|
|
|
|
|
def _log_map(t, c, s):
|
|
"""u = sign(t - c) * log1p(|t - c| / s). An ADDITIVE step in u is a
|
|
MULTIPLICATIVE step in x, which is what the wide-bracket families need: the
|
|
search there must cross up to twelve decades, and an absolute step sized by
|
|
the local slope knows nothing about how many decades remain."""
|
|
d = t - c
|
|
return math.copysign(math.log1p(abs(d) / s), d)
|
|
|
|
|
|
def _log_unmap(u, c, s):
|
|
return c + math.copysign(s * (math.expm1(abs(u))), u)
|
|
|
|
|
|
def _log_space_step(S, x, fx, far, slope, c, s):
|
|
"""The same similarity rule, the same derived scale, read in u rather than x.
|
|
|
|
dx/du = s + |x - c|, so the slope in u is slope_x * (s + |x - c|); the scale
|
|
is that slope times the distance to the far end IN u, exactly as in x-space.
|
|
Returns the step in x (None when u-space cannot be read).
|
|
|
|
`s` is NOT a free unit: PREREG amendment 1 fixes s = xtol by derivation - the
|
|
only length in the problem that says "below this nothing is distinguishable",
|
|
supplied by the caller like a, b and f. With it u spans log(w/xtol) units,
|
|
the very measure the method's own guarantee 2*ceil(log2(w/xtol)) + 3 is
|
|
stated in, so the march is linear in the measure of the guarantee.
|
|
"""
|
|
u_x = _log_map(x, c, s)
|
|
u_far = _log_map(far, c, s)
|
|
du = u_far - u_x
|
|
if du == 0.0 or not math.isfinite(du):
|
|
return None
|
|
slope_u = slope * (s + abs(x - c))
|
|
sc_u = abs(slope_u) * abs(du)
|
|
if sc_u == 0.0 or not math.isfinite(sc_u):
|
|
return None
|
|
u_new = u_x + (1.0 - S(abs(fx) / sc_u)) * du
|
|
x_new = _log_unmap(u_new, c, s)
|
|
if not math.isfinite(x_new):
|
|
return None
|
|
return x_new - x
|
|
|
|
|
|
class _Safeguard:
|
|
"""Parameter-free worst-case protection. See the module docstring."""
|
|
|
|
__slots__ = ("mode", "guard_width", "pending", "step_2", "step_1",
|
|
"n_half", "eps", "j", "allow", "credit")
|
|
|
|
def __init__(self, mode, width, xtol=None, budget_factor=2.0,
|
|
budget_extra=None):
|
|
if mode not in GUARDS:
|
|
raise ValueError(f"unknown guard {mode!r}, expected one of {GUARDS}")
|
|
self.mode = mode
|
|
self.guard_width = width
|
|
self.pending = False
|
|
self.step_2 = width # step two evaluations ago
|
|
self.step_1 = width # step one evaluation ago
|
|
# projection state: the bisection budget, fixed by the initial bracket
|
|
self.eps = 0.5 * xtol if xtol else 0.0
|
|
self.j = 0
|
|
# set for every mode: the step loop decrements it unconditionally, and
|
|
# leaving it to the projection branch made every other guard raise
|
|
self.allow = 0.0
|
|
self.n_half = 0
|
|
self.credit = 1.0 # one bit of endowment; see project_credit
|
|
if mode in ("project", "safe"):
|
|
if not xtol or width <= 0.0:
|
|
raise ValueError("the projection guard needs xtol and a "
|
|
"positive initial width")
|
|
half = math.ceil(math.log2(width / xtol)) if width > xtol else 0
|
|
# The budget is TWICE the bisection budget. Zero slack (budget =
|
|
# half) is degenerate: a worst case of exactly bisection demands a
|
|
# halving at every single step, so the projection leaves no
|
|
# freedom at all - measured, r_0 was 0.05 of a unit bracket and
|
|
# the guard behaved as bisection. Doubling is not a tuned
|
|
# constant: it is the same factor 2 already claimed as the bound,
|
|
# and it buys total freedom until the budget is nearly spent.
|
|
# Two ways to give the projection room, and the difference
|
|
# matters for whether this method contains a constant at all.
|
|
#
|
|
# budget_extra = k -> n_max = half + k. The requirement is
|
|
# only r_0 >= 0, i.e. n_max >= half; at n_max = half there
|
|
# is no room and the guard degenerates to bisection. The
|
|
# MINIMAL non-degenerate choice is therefore k = 1 - not a
|
|
# tuned value but the smallest integer for which the
|
|
# construction works at all.
|
|
# budget_factor = c -> n_max = ceil(c * half). This is what
|
|
# the paper used, at c = 2, and the justification given for
|
|
# it ("the factor already in the bound") was circular: the
|
|
# bound is 2x because 2 was chosen. Kept as an ablation
|
|
# knob so the circularity can be measured rather than
|
|
# argued about.
|
|
if budget_extra is not None:
|
|
self.n_half = half + int(budget_extra)
|
|
else:
|
|
self.n_half = max(half, int(math.ceil(budget_factor * half)))
|
|
# The allowance 2^(n_max - j).eps is halved once per step, so it
|
|
# is carried forward rather than recomputed: an exponentiation in
|
|
# the inner loop was costing more than the projection itself.
|
|
# Halving is exact in binary floating point, so this changes no
|
|
# iterate and no evaluation count - only the arithmetic.
|
|
try:
|
|
self.allow = 2.0 ** self.n_half * self.eps
|
|
except OverflowError:
|
|
self.allow = math.inf
|
|
|
|
def project_credit(self, x_sim, lo, hi):
|
|
"""The projection with the budget EARNED rather than scheduled.
|
|
|
|
The scheduled version fixes the allowance in advance from the initial
|
|
width, so the multiplier in n_max = c*half is a free constant, and the
|
|
argument that c = 2 is "the factor already in the bound" is circular:
|
|
the bound is 2x because 2 was chosen. Measured, c = 1.25 is 8%
|
|
better, so the constant is not neutral either.
|
|
|
|
Here the method pays its own way. Let s_j = log2(w_j / w_{j+1}) be
|
|
the shrink a step delivers in bits; bisection delivers exactly 1. Keep
|
|
|
|
C <- C + s_j - 1, C >= 0,
|
|
|
|
the surplus banked so far, and allow a step to under-deliver only by
|
|
what has been banked: the width after the step may be at most
|
|
2^(C-1) * w. Since C never goes negative, the shrinks over j steps
|
|
sum to at least j - C_0, so the run costs at most half + C_0 steps.
|
|
|
|
C_0 = 1 - one bit of endowment - is the smallest integer for which the
|
|
rule is not degenerate: at C_0 = 0 the allowance is exactly w/2 on the
|
|
first step, which forces the midpoint and makes the method bisection.
|
|
So the guarantee is bisection plus one step, and the constant is the
|
|
promise, not a tuning knob.
|
|
"""
|
|
w = hi - lo
|
|
# The allowance is w * 2^(C-1), so it exceeds w exactly when C >= 1,
|
|
# and then it restricts nothing. Testing that FIRST removes the
|
|
# exponentiation from the common path and, more importantly, removes
|
|
# a crash: C accumulates the log of every shrink and is unbounded, so
|
|
# 2**(C-1) overflows a double at C ~ 1025. That is reachable - sixty
|
|
# nine-decade steps reach C = 1012, and a bracket of width 1e300 with
|
|
# a tolerance of 1e-8 has about 1023 bits of shrink available - and
|
|
# the OverflowError would escape, since only the budget exception is
|
|
# caught here.
|
|
if self.credit >= 1.0:
|
|
return x_sim
|
|
allowed = w * (2.0 ** (self.credit - 1.0))
|
|
left, right = hi - allowed, lo + allowed
|
|
if left > right:
|
|
return 0.5 * (lo + hi)
|
|
return min(max(x_sim, left), right)
|
|
|
|
def note_shrink(self, w_before, w_after):
|
|
"""Bank what the step delivered beyond a halving, or draw on the bank
|
|
for what it fell short."""
|
|
if w_before <= 0.0 or w_after <= 0.0:
|
|
return
|
|
s = math.log2(w_before / w_after)
|
|
self.credit = max(0.0, self.credit + s - 1.0)
|
|
|
|
def project(self, x_sim, lo, hi):
|
|
"""Pull the similarity step back to the edge of the interval that the
|
|
bisection budget still allows. Returns the admissible iterate.
|
|
|
|
The budget says the next bracket may be at most
|
|
`allowed = 2^(n_max - j) * eps`. Whichever side the sign change turns
|
|
out to be on, that bracket is at most max(x - lo, hi - x), so the
|
|
admissible set is exactly
|
|
|
|
[hi - allowed, lo + allowed].
|
|
|
|
Writing it this way rather than as a radius r = allowed - (hi-lo)/2
|
|
around the midpoint matters numerically: r is the difference of two
|
|
nearly equal numbers, and the cancellation made the computed radius
|
|
slightly too large. Measured, that admitted points which broke the
|
|
invariant by a relative 8.9e-5 and cost exactly one extra evaluation
|
|
on eleven Chandrupatla problems. The clamped form below performs no
|
|
cancelling subtraction and cannot overshoot."""
|
|
allowed = self.allow
|
|
if not math.isfinite(allowed):
|
|
return x_sim
|
|
left, right = hi - allowed, lo + allowed
|
|
if left > right:
|
|
# the budget no longer permits anything but a halving
|
|
return 0.5 * (lo + hi)
|
|
return min(max(x_sim, left), right)
|
|
|
|
def bisect_before_step(self, width):
|
|
"""Called before the similarity step is even computed. Only the
|
|
bracket guard can decide at this point."""
|
|
if self.mode not in ("bracket", "both"):
|
|
return False
|
|
if width <= 0.5 * self.guard_width:
|
|
self.guard_width = width
|
|
# The similarity step that is about to be taken is itself the
|
|
# pending one. Clearing the flag here instead granted TWO
|
|
# consecutive similarity steps after every halving, which makes
|
|
# the guarantee one bisection in three evaluations (3x bisection)
|
|
# rather than one in two (2x) - measured at 139 evaluations
|
|
# against a 2x bound of 114 on Chandrupatla fun4.3.
|
|
self.pending = True
|
|
return False
|
|
if self.pending:
|
|
self.pending = False
|
|
return True
|
|
self.pending = True
|
|
return False
|
|
|
|
def reject_step(self, step):
|
|
"""Called with the proposed displacement. Only the step guard can
|
|
decide at this point."""
|
|
if self.mode not in ("step", "both", "safe"):
|
|
return False
|
|
return abs(step) > 0.5 * self.step_2
|
|
|
|
def note_step(self, step):
|
|
self.step_2 = self.step_1
|
|
self.step_1 = abs(step)
|
|
|
|
|
|
def _slope3(x2, f2, x1, f1, x0, f0):
|
|
"""Derivative at x2 of the parabola through the three most recent points,
|
|
as Newton divided differences:
|
|
|
|
f'(x2) ~ f[x1,x2] + f[x0,x1,x2] * (x2 - x1)
|
|
|
|
Falls back to the two-point slope when a point is repeated or missing, so
|
|
the first iteration of a run behaves exactly like scale="derived"."""
|
|
if x2 == x1:
|
|
return 0.0
|
|
d12 = (f2 - f1) / (x2 - x1)
|
|
# An overflowing difference must report "no usable slope" (0.0), which
|
|
# sends pure_root to its bisection fallback. Returning the infinity
|
|
# instead would make the scale infinite, hence z = 0, hence "perfectly
|
|
# similar to the target", hence a zero step - a stall dressed up as
|
|
# success.
|
|
if not math.isfinite(d12):
|
|
return 0.0
|
|
if x0 is None or x0 == x1 or x0 == x2:
|
|
return d12
|
|
d01 = (f1 - f0) / (x1 - x0)
|
|
d012 = (d12 - d01) / (x2 - x0)
|
|
s = d12 + d012 * (x2 - x1)
|
|
return s if math.isfinite(s) else d12
|
|
|
|
|
|
def _local_order(x2, f2, x1, f1, x0, f0):
|
|
"""The local order m of f, from three evaluations and nothing else.
|
|
|
|
Redo the contraction law WITHOUT assuming |f| ~ |f'(r)|.|x-r|. Take
|
|
instead |f| ~ K.e^m near the current point, m = 1 being a simple root and
|
|
m > 1 a locally flat approach. Then for a power similarity S = 1 - z^p,
|
|
|
|
z = |f| / (|f'|.D) = e / (m.D) (K cancels)
|
|
step = z^p . D = D^(1-p) . (e/m)^p
|
|
|
|
and the step lands on the root when (e/D)^(p-1) = m^p, i.e.
|
|
|
|
p = L / (L - ln m), L = ln(e/D) < 0.
|
|
|
|
m = 1 gives p = 1, the earlier law. m > 1 drives p below 1 on its own -
|
|
at e/D = 0.1, m = 10 it gives p = 0.50, and at e/D = 0.5 it gives 0.23,
|
|
which is the neighbourhood of the cube-root similarity that stage 7 found
|
|
winning on the wide-bracket families. Substituting that p back collapses
|
|
the whole expression to step = e = m.|f|/|f'|, so the exponent is a
|
|
reparameterisation and the substance is the multiplier m.
|
|
|
|
m is measurable without knowing the root: with theta = f.f''/(f')^2,
|
|
f ~ K e^m, f' ~ K m e^(m-1), f'' ~ K m(m-1) e^(m-2)
|
|
=> theta = m(m-1)/m^2 = 1 - 1/m => m = 1/(1 - theta).
|
|
Both derivatives are divided differences of the three most recent points,
|
|
so nothing is chosen or fitted.
|
|
|
|
Returns 1.0 whenever the power-law model does not apply - theta outside
|
|
(0,1) means the data does not describe a flat approach, and the honest
|
|
response is to decline the correction rather than to clamp it into range.
|
|
"""
|
|
if x0 is None or len({x0, x1, x2}) < 3:
|
|
return 1.0
|
|
d01 = (f1 - f0) / (x1 - x0)
|
|
d12 = (f2 - f1) / (x2 - x1)
|
|
d012 = (d12 - d01) / (x2 - x0)
|
|
fp = d12 + d012 * (x2 - x1) # f'(x2) on the parabola
|
|
fpp = 2.0 * d012 # f''(x2) on the parabola
|
|
if fp == 0.0 or not math.isfinite(fp) or not math.isfinite(fpp):
|
|
return 1.0
|
|
theta = f2 * fpp / (fp * fp)
|
|
if not math.isfinite(theta) or not (0.0 < theta < 1.0):
|
|
return 1.0
|
|
m = 1.0 / (1.0 - theta)
|
|
return m if math.isfinite(m) and m >= 1.0 else 1.0
|
|
|
|
|
|
def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived",
|
|
guard="step", stats=None, budget_factor=2.0,
|
|
budget_extra=None, endowment=1.0, space="x", prestep="none",
|
|
points="recent", direction="distance", exponent="none"):
|
|
"""The pure similarity root finder.
|
|
|
|
similarity : key into PURE_SIMS
|
|
scale : "derived" - local slope of f from the two most recent
|
|
evaluations (parameter free, order ~1.618)
|
|
"derived3_m" - REFUTED as an improvement. The three-point
|
|
scale divided by the local order m of f, which
|
|
is what the contraction law gives once the
|
|
linearity assumption is dropped. The derivation
|
|
is sound and the order is preserved; the
|
|
estimator of m is not sharp enough to pay for
|
|
itself. See the note at the use site.
|
|
"derived_rho", "derived3_rho" - REFUTED. The same, with the
|
|
scale multiplied by the fraction of its own
|
|
promise the previous step kept. Retained for
|
|
reproducibility only; see the note at the use
|
|
site. Do not adopt.
|
|
"derived3" - local slope from the three most recent
|
|
evaluations, i.e. the derivative at x of the
|
|
parabola through them. Still parameter free -
|
|
a divided difference is a computation, not a
|
|
fitted coefficient - but the slope error drops
|
|
from O(h) to O(h^2), which lifts the order from
|
|
1.618 towards 1.84. This is the direct
|
|
consequence of the stage-4 finding that the
|
|
whole deficit against Brent is the endgame
|
|
order, and nothing else.
|
|
"start" - |f| at the starting point (a chosen reference)
|
|
"bracket" - largest |f| at the current bracket ends
|
|
guard : one of GUARDS
|
|
endowment : ABLATION KNOB for guard="credit", not a setting. The method
|
|
is defined at 1 bit, and 1 is not a chosen value: the
|
|
guarantee is bisection + ceil(C0) steps, so every C0 in (0,1]
|
|
promises the same thing - one step - and among those 1 is the
|
|
most permissive and therefore the cheapest. Measured over
|
|
the 199 problems: 0.1 -> 8395, 0.5 -> 7657, 1 -> 6588, all
|
|
with the same one-step guarantee; 2 -> 6404 but two steps,
|
|
5 -> 5621 but five. Nothing is free to choose once the
|
|
promise is fixed, and the promise is the specification.
|
|
stats : optional dict; on return it carries 'forced' = how many steps
|
|
the safeguard turned into bisections. Comparable with the
|
|
same key from baselines.brent, so the two methods' fallback
|
|
rates can be put side by side.
|
|
|
|
Returns (x, nfev, converged).
|
|
"""
|
|
if similarity not in PURE_SIMS:
|
|
raise ValueError(f"unknown similarity {similarity!r}")
|
|
if scale not in ("derived", "derived3", "derived_rho", "derived3_rho",
|
|
"derived3_m", "derived3_cap", "start", "bracket", "far"):
|
|
raise ValueError(f"unknown scale {scale!r}")
|
|
if exponent not in EXPONENTS:
|
|
raise ValueError(f"unknown exponent {exponent!r}")
|
|
if exponent == "order" and scale not in ("derived3", "derived3_cap", "bracket", "far"):
|
|
# the order is read from a three-point divided difference; under the two-point
|
|
# scale there is no second difference to read it from, and computing one from
|
|
# elsewhere would hide that (PREREG 2026-09-21, slot 1)
|
|
raise ValueError("exponent='order' requires scale in "
|
|
"{'derived3', 'derived3_cap', 'bracket', 'far'}")
|
|
if space not in SPACES:
|
|
raise ValueError(f"unknown space {space!r}")
|
|
if prestep not in ("none", "progress", "slope", "band", "defect", "simpair"):
|
|
raise ValueError(f"unknown prestep {prestep!r}")
|
|
if prestep in ("band", "defect", "simpair") and scale not in (
|
|
"derived3", "derived3_cap", "derived3_m", "bracket", "far"):
|
|
# the decision reads a TRIPLE; under the two-point scale there is no third point
|
|
# to read a configuration from, and reading one from elsewhere would hide that
|
|
raise ValueError("prestep in {'band','defect','simpair'} needs a three-point "
|
|
"history (scale derived3*/bracket/far)")
|
|
if points not in POINT_RULES:
|
|
raise ValueError(f"unknown points rule {points!r}")
|
|
if direction not in DIRECTIONS:
|
|
raise ValueError(f"unknown direction {direction!r}")
|
|
if points == "shape" and scale not in ("derived3", "derived3_rho", "derived3_m",
|
|
"derived3_cap"):
|
|
# the triple only exists for the three-point scale; selecting one under a
|
|
# two-point scale would change nothing and hide that
|
|
raise ValueError("points='shape' requires a three-point scale (derived3*)")
|
|
if prestep == "slope" and scale not in ("derived", "derived3", "derived3_cap"):
|
|
# there is no slope to judge under a chosen or bracket scale; judging something
|
|
# else would hide that, so the call fails instead
|
|
raise ValueError("prestep='slope' requires scale in {'derived', 'derived3'}")
|
|
if space == "log" and scale not in ("derived", "derived3"):
|
|
# u-space needs a slope; with a chosen or bracket scale there is none to
|
|
# transform, and computing something anyway would hide that (PREREG §3a).
|
|
raise ValueError("space='log' requires scale in {'derived', 'derived3'}")
|
|
S = PURE_SIMS[similarity]
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b)
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
if fa * fb > 0.0:
|
|
raise ValueError("f(a) and f(b) must straddle zero")
|
|
|
|
A, FA, B, FB = a, fa, b, fb # the proven bracket
|
|
if abs(fa) <= abs(fb):
|
|
x, fx, x_prev, f_prev = a, fa, b, fb
|
|
else:
|
|
x, fx, x_prev, f_prev = b, fb, a, fa
|
|
# the third point does not exist yet; _slope3 falls back to two points
|
|
x_prev2, f_prev2 = None, None
|
|
start_scale = abs(fx)
|
|
sg = _Safeguard(guard, abs(B - A), xtol, budget_factor, budget_extra)
|
|
sg.credit = endowment
|
|
forced = 0
|
|
pre_refused = 0
|
|
triple_changed = 0
|
|
dir_changed = 0
|
|
slope = 0.0
|
|
slope_prev = None
|
|
if stats is not None:
|
|
stats['forced'] = 0
|
|
stats['pre_refused'] = 0
|
|
stats['triple_changed'] = 0
|
|
stats['dir_changed'] = 0
|
|
|
|
while abs(B - A) > xtol:
|
|
if sg.bisect_before_step(abs(B - A)):
|
|
x_new = 0.5 * (A + B)
|
|
forced += 1
|
|
# keep the step history honest under guard="both": a bisection
|
|
# forced by the bracket backstop is still a step that was taken
|
|
sg.note_step(x_new - x)
|
|
else:
|
|
far = B if abs(B - x) >= abs(A - x) else A
|
|
if direction == "dissimilar":
|
|
# PREREG 2026-09-21, mode B: the far end read by SIMILARITY instead of
|
|
# geometry — the end whose |f| is larger is the one LESS similar to zero,
|
|
# i.e. the one the root is further from. Two values already in hand are
|
|
# compared; no constant.
|
|
far_sim = B if abs(FB) >= abs(FA) else A
|
|
if far_sim != far:
|
|
dir_changed += 1
|
|
far = far_sim
|
|
if scale in ("derived", "derived_rho"):
|
|
slope = ((fx - f_prev) / (x - x_prev)) if x != x_prev else 0.0
|
|
sc = abs(slope) * abs(far - x)
|
|
elif scale in ("derived3", "derived3_rho", "derived3_m"):
|
|
if points == "shape":
|
|
# PREREG 2026-09-21, mode A: which THREE points the model is built
|
|
# on, chosen by the SHAPE of the configuration instead of by recency.
|
|
# For a linear f the position ratio and the value ratio coincide, so
|
|
# |phi - xi| is a dimensionless measure of how unlike a linear
|
|
# configuration a triple is; the triple minimising it is SELECTED
|
|
# (the corpus rule is select, never weight). argmin over the declared
|
|
# candidates needs no threshold.
|
|
pick = _pick_triple(x, fx, x_prev, f_prev, x_prev2, f_prev2,
|
|
A, FA, B, FB, far)
|
|
if pick is None:
|
|
slope = _slope3(x, fx, x_prev, f_prev, x_prev2, f_prev2)
|
|
else:
|
|
(xu, fu), (xv, fv) = pick
|
|
if xu != x_prev or xv != x_prev2:
|
|
triple_changed += 1
|
|
slope = _slope3(x, fx, xu, fu, xv, fv)
|
|
else:
|
|
slope = _slope3(x, fx, x_prev, f_prev, x_prev2, f_prev2)
|
|
sc = abs(slope) * abs(far - x)
|
|
if scale == "derived3_m":
|
|
# REFUTED as an improvement; kept so the refutation stays
|
|
# REPRODUCIBLE, which is the whole point of keeping it. The long
|
|
# note is at the end of the derived3_cap branch below.
|
|
sc /= _local_order(x, fx, x_prev, f_prev, x_prev2, f_prev2)
|
|
elif scale == "derived3_cap":
|
|
# Valeri, 2026-09-21: "if the similarity formed the step, the step would be
|
|
# large far from the target and small near it". It does not, and the reason
|
|
# is the normalisation: with sigma = |slope| * D the similarity reads "what
|
|
# fraction of the BRACKET the Newton step is", not "how far am I". On a
|
|
# 2e4-wide bracket of x^3-2x-5 at x = 1e4 that fraction is 0.167, so S =
|
|
# 0.833 - "similar" - and the step is small exactly where it should be large.
|
|
# The scale is therefore SELECTED between the two already-computed ones:
|
|
# far away the value scale of the bracket ends is the smaller and the step
|
|
# opens up; near the root the slope scale is the smaller and Newton (and the
|
|
# order) returns. A min over what is in hand - no constant, no threshold.
|
|
slope = _slope3(x, fx, x_prev, f_prev, x_prev2, f_prev2)
|
|
sc = min(abs(slope) * abs(far - x), max(abs(FA), abs(FB)))
|
|
# NOTE on scale="derived3_m" (the live code is in the derived3 branch
|
|
# above). On 2026-09-21 the derived3_cap branch was inserted directly
|
|
# in front of this note and swallowed the division with it, so
|
|
# derived3_m silently became derived3 - 2806 against 2806 - and the
|
|
# refutation stopped being reproducible. test_m_loses_on_the_full
|
|
# _benchmark caught it, which is what that pin is for.
|
|
# REFUTED as an improvement - kept for reproducibility.
|
|
# The DERIVATION holds: dropping the linearity
|
|
# assumption from the contraction law does produce the
|
|
# march exponent, and it collapses to the multiplier m.
|
|
# The ESTIMATOR does not: theta from three points
|
|
# cannot tell a genuinely flat approach from curvature
|
|
# noise. Measured, it reports m > 1 on 27.9% of steps
|
|
# across Alefeld-Potra-Shi, where almost every root is
|
|
# simple and the truth is m = 1, with excursions to
|
|
# m = 34.8; on the flat fun3/4/5 families it fires on
|
|
# 35.8%. A separation of 27.9 against 35.8 is no
|
|
# separation, and the false positives cost 273
|
|
# evaluations: 3079 against 2806 on APS, 6304 against
|
|
# 6026 over both sets.
|
|
#
|
|
# It does pass the two tests rho failed - the order
|
|
# stays at 2.001 and the eight standard problems
|
|
# improve to 82 - which is exactly why it must be
|
|
# labelled: a small hand-picked set makes it look like
|
|
# a win.
|
|
#
|
|
# step = m.|f|/|f'| is step = e for a locally flat
|
|
# approach; dividing the scale by m is the same thing
|
|
# written in the similarity's own terms. m >= 1, so
|
|
# like rho this can only lengthen a step - but unlike
|
|
# rho it is read from the CURVATURE, which stays
|
|
# meaningful as |f| goes to zero, and it returns to 1
|
|
# at a simple root, leaving the endgame untouched.
|
|
elif scale == "start":
|
|
sc = start_scale
|
|
elif scale == "far":
|
|
# the value at the end the step AIMS AT, already computed. With
|
|
# S = 1/(1+z) this is false position exactly (Proposition 4), which
|
|
# is why it is the honest control for the order exponent: the pair
|
|
# (far, none) against (far, order) differs in one thing only.
|
|
sc = abs(FB) if far == B else abs(FA)
|
|
else:
|
|
sc = max(abs(FA), abs(FB))
|
|
if scale.endswith("_rho"):
|
|
# REFUTED - kept only so the experiment is reproducible.
|
|
# Measured (stage 8): the fun3/4/5 families it was built
|
|
# for do not improve (1364 against 1367, where cbrt gets
|
|
# 677); the eight standard problems get twice as
|
|
# expensive (84 -> 161); and the observed order collapses
|
|
# from 2.000 to 0.904.
|
|
#
|
|
# Why the derivation was wrong: rho equates "how much
|
|
# value the step removed" with "how trustworthy the local
|
|
# model is". Near the root |f_n| and |f_{n-1}| are both
|
|
# negligible and their ratio is noise, so rho thrashes in
|
|
# exactly the place where the correction has to be inert.
|
|
# The alternative credibility - the ratio of the assumed
|
|
# slope to the realised one - is closed too: the realised
|
|
# slope is what the derived scale already computes, so
|
|
# correcting by it would double-count.
|
|
#
|
|
# The phase detector. A step of |f|/|f'| is a promise to
|
|
# land on the root; rho is the fraction of that promise
|
|
# the previous step actually kept:
|
|
#
|
|
# rho = 1 - |f_n| / |f_{n-1}|
|
|
#
|
|
# Near the root the local model is right and rho -> 1, so
|
|
# the correction vanishes and the order is untouched. Far
|
|
# from it - the approach phase on a very wide bracket -
|
|
# the model under-delivers, rho is small, and the step
|
|
# must be 1/rho longer to remove what is left. Dividing
|
|
# the scale by rho does exactly that.
|
|
#
|
|
# rho <= 1 always, so the correction can only LENGTHEN the
|
|
# step, never shorten it below the Newton step. rho <= 0
|
|
# means the value did not improve at all: there is no
|
|
# credibility to read, so the correction is dropped rather
|
|
# than guessed, and the bracket projection keeps the step
|
|
# safe either way.
|
|
if f_prev != 0.0:
|
|
rho = 1.0 - abs(fx) / abs(f_prev)
|
|
if 0.0 < rho <= 1.0:
|
|
sc *= rho
|
|
if sc == 0.0 or not math.isfinite(sc):
|
|
x_new = 0.5 * (A + B)
|
|
else:
|
|
m_loc = (_local_order(x, fx, x_prev, f_prev, x_prev2, f_prev2)
|
|
if exponent == "order" else None)
|
|
z_now = _read_z(abs(fx), sc, exponent, m_loc)
|
|
step = None
|
|
if prestep in ("band", "defect", "simpair"):
|
|
# THE DECISION PLACE (PREREG 2026-09-21, slot 2). Same slot as the
|
|
# slope gate, a different reading of it: the question is not whether
|
|
# the slope has settled but whether the CONFIGURATION deserves the
|
|
# model. The triple is the recency triple, so only the reading
|
|
# changes against the slope-gate cells.
|
|
if _decide_refuse(prestep, S, z_now, x, fx, x_prev, f_prev,
|
|
x_prev2, f_prev2):
|
|
step = 0.5 * (A + B) - x
|
|
pre_refused += 1
|
|
if step is None and prestep == "slope":
|
|
# PRE-STEP GATE ON SHAPE (PREREG 2026-09-21). The step trusts the
|
|
# local slope, so the question of trust is whether that slope has
|
|
# SETTLED: two consecutive estimates of opposite sign mean the local
|
|
# model is not established at all, and the evaluation is better spent
|
|
# halving. A SIGN is compared, not a magnitude — no threshold, no
|
|
# factor, no constant. Where the slope sign is stable (a simple root
|
|
# of a smooth f) the gate is inert, which is exactly what the
|
|
# progress gate could not manage: it charged APS 4449 against 2806.
|
|
if (slope_prev is not None and slope != 0.0 and slope_prev != 0.0
|
|
and (slope > 0.0) != (slope_prev > 0.0)):
|
|
step = 0.5 * (A + B) - x
|
|
pre_refused += 1
|
|
if step is None and prestep == "progress":
|
|
# PRE-STEP GATE (PREREG 2026-09-21). Measured cause of our cost:
|
|
# the safeguard REPAIRS 77.5% of evaluations after the fact, while
|
|
# Chandrupatla REFUSES its model on 72% of steps for free and
|
|
# halves instead - 2.310 halvings per evaluation against our 1.422.
|
|
# Our step never leaves the bracket (0 <= S <= 1), so the question
|
|
# is not overshoot but PROGRESS: if the similarity step promises
|
|
# less than bisection would remove, spend the same single
|
|
# evaluation on the midpoint instead. Two lengths already in hand
|
|
# are compared; there is no threshold and no new constant.
|
|
cand = (1.0 - S(z_now)) * (far - x)
|
|
half = 0.5 * (A + B) - x
|
|
if abs(cand) < abs(half):
|
|
step = half
|
|
pre_refused += 1
|
|
if step is None and space == "log":
|
|
# c = the bracket end with the smaller |f| (the side the
|
|
# root is nearer); s = xtol, fixed by derivation in
|
|
# PREREG amendment 1, never chosen. Falls back to the
|
|
# x-space step when u cannot be read, so the mode can
|
|
# only differ where it is defined.
|
|
c_near = A if abs(FA) <= abs(FB) else B
|
|
step = _log_space_step(S, x, fx, far, slope, c_near, xtol)
|
|
if step is None:
|
|
step = (1.0 - S(z_now)) * (far - x)
|
|
if sg.reject_step(step):
|
|
x_new = 0.5 * (A + B)
|
|
step = x_new - x
|
|
forced += 1
|
|
else:
|
|
x_new = x + step
|
|
if guard in ("project", "safe", "credit"):
|
|
projected = (sg.project_credit if guard == "credit"
|
|
else sg.project)(x_new, min(A, B), max(A, B))
|
|
if projected != x_new:
|
|
forced += 1
|
|
x_new = projected
|
|
step = x_new - x
|
|
sg.note_step(step)
|
|
# the slope this step trusted, kept so the NEXT step can ask whether it
|
|
# has settled (prestep="slope"); only meaningful for the derived scales
|
|
slope_prev = slope
|
|
|
|
if not math.isfinite(x_new) or not (min(A, B) < x_new < max(A, B)):
|
|
x_new = 0.5 * (A + B)
|
|
if x_new == x or x_new == x_prev:
|
|
x_new = 0.5 * (A + B)
|
|
if x_new == x or x_new == x_prev:
|
|
# the bracket can no longer be split in floating point
|
|
if stats is not None:
|
|
stats['forced'] = forced; stats['pre_refused'] = pre_refused; stats['triple_changed'] = triple_changed; stats['dir_changed'] = dir_changed
|
|
return x_new, g.n, abs(B - A) <= xtol
|
|
|
|
w_before = abs(B - A)
|
|
fnew = g(x_new)
|
|
sg.j += 1 # one step of the bisection budget is spent
|
|
sg.allow *= 0.5
|
|
x_prev2, f_prev2 = x_prev, f_prev
|
|
x_prev, f_prev = x, fx
|
|
x, fx = x_new, fnew
|
|
if fnew == 0.0:
|
|
if stats is not None:
|
|
stats['forced'] = forced; stats['pre_refused'] = pre_refused; stats['triple_changed'] = triple_changed; stats['dir_changed'] = dir_changed
|
|
return x_new, g.n, True
|
|
if FA * fnew < 0.0:
|
|
B, FB = x_new, fnew
|
|
else:
|
|
A, FA = x_new, fnew
|
|
if guard == "credit":
|
|
sg.note_shrink(w_before, abs(B - A))
|
|
if stats is not None:
|
|
stats['forced'] = forced; stats['pre_refused'] = pre_refused; stats['triple_changed'] = triple_changed; stats['dir_changed'] = dir_changed
|
|
return 0.5 * (A + B), g.n, True
|
|
except _Budget:
|
|
if stats is not None:
|
|
stats['forced'] = forced; stats['pre_refused'] = pre_refused; stats['triple_changed'] = triple_changed; stats['dir_changed'] = dir_changed
|
|
return x, g.n, False
|