695 lines
35 KiB
Python
695 lines
35 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")
|
|
|
|
|
|
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"):
|
|
"""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", "start", "bracket"):
|
|
raise ValueError(f"unknown scale {scale!r}")
|
|
if space not in SPACES:
|
|
raise ValueError(f"unknown space {space!r}")
|
|
if prestep not in ("none", "progress", "slope"):
|
|
raise ValueError(f"unknown prestep {prestep!r}")
|
|
if prestep == "slope" and scale not in ("derived", "derived3"):
|
|
# 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
|
|
slope = 0.0
|
|
slope_prev = None
|
|
if stats is not None:
|
|
stats['forced'] = 0
|
|
stats['pre_refused'] = 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 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"):
|
|
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 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.
|
|
sc /= _local_order(x, fx, x_prev, f_prev,
|
|
x_prev2, f_prev2)
|
|
elif scale == "start":
|
|
sc = start_scale
|
|
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:
|
|
step = None
|
|
if 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(abs(fx) / sc)) * (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(abs(fx) / sc)) * (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
|
|
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
|
|
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
|
|
return 0.5 * (A + B), g.n, True
|
|
except _Budget:
|
|
if stats is not None:
|
|
stats['forced'] = forced; stats['pre_refused'] = pre_refused
|
|
return x, g.n, False
|