# cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False # cython: embedsignature=False, binding=False, always_allow_keywords=True # # The compiled A' bracketing root finder. This file is the SOURCE of the distributed library # and is not shipped; the wheel carries only the compiled extension. # # Port of simroot.pure.pure_root restricted to the one configuration the paper proposes: # similarity="ratio", scale="derived3", guard="project" (budget factor 2), space="x", # prestep="contra+flatanalogmempt+exp+half+close+sub", points="recent", # direction="distance", exponent="none", extrapolate="pattern". # It must reproduce pure_root bit for bit (x and the evaluation count) - pinned by # tests/test_simbracket_equivalence.py on every problem of the six corpora. Python's own # semantics are kept wherever they differ from C: min/max keep the FIRST argument on ties and # NaN, float division by zero raises ZeroDivisionError (cdivision stays off), math.log raises # ValueError at 0, and the two library calls whose result depends on the interpreter # (math.hypot and the builtin sum, compensated since Python 3.12) are called, not re-coded. from libc.math cimport fabs, log, log2, exp, pow, sqrt, copysign, isfinite, isnan, isinf from libc.math cimport ceil, INFINITY from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free import builtins as _builtins from math import hypot as _py_hypot _py_sum = _builtins.sum __all__ = ["find_root", "evaluation_bound"] class _Budget(Exception): pass # ---------------------------------------------------------------------------------------- # Python semantics of the scalar helpers # ---------------------------------------------------------------------------------------- cdef inline double _pmax(double a, double b) noexcept: # builtins.max(a, b): the first argument unless the second is strictly greater return b if b > a else a cdef inline double _pmin(double a, double b) noexcept: # builtins.min(a, b): the first argument unless the second is strictly smaller return b if b < a else a cdef inline double _sim(double z) noexcept: return _pmax(0.0, 1.0 - fabs(z)) cdef inline bint _opposite(double u, double v) noexcept: return (u < 0.0 < v) or (v < 0.0 < u) cdef inline bint _same_sign(double u, double v) noexcept: return (u > 0.0 and v > 0.0) or (u < 0.0 and v < 0.0) cdef double _plog(double v) except? -1.5: # math.log: nan -> nan, +inf -> +inf, v <= 0 or -inf -> ValueError if v > 0.0 or isnan(v): return log(v) raise ValueError("math domain error") cdef inline double _ratio_sim(double a, double b) except? -1.5: cdef double hi = _pmax(a, b) return _pmin(a, b) / hi if hi > 0.0 else 0.0 # ---------------------------------------------------------------------------------------- # growable double buffer # ---------------------------------------------------------------------------------------- cdef struct Buf: double* p Py_ssize_t n Py_ssize_t cap cdef int _push(Buf* b, double v) except -1: cdef Py_ssize_t cap cdef double* q if b.n == b.cap: cap = 64 if b.cap == 0 else 2 * b.cap q = PyMem_Realloc(b.p, cap * sizeof(double)) if q == NULL: raise MemoryError() b.p = q b.cap = cap b.p[b.n] = v b.n += 1 return 0 # ---------------------------------------------------------------------------------------- # the pieces of the step # ---------------------------------------------------------------------------------------- cdef double _slope3(double x2, double f2, double x1, double f1, bint has0, double x0, double f0) except? -1.5: cdef double d12, d01, d012, s if x2 == x1: return 0.0 d12 = (f2 - f1) / (x2 - x1) if not isfinite(d12): return 0.0 if not has0 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 isfinite(s) else d12 cdef double _order_split(double A, double FA, double B, double FB, double m) except? -1.5: cdef double lq, t, e if isinf(m): return 0.5 * (A + B) lq = (_plog(fabs(FA)) - _plog(fabs(FB))) / m if not isfinite(lq): return 0.5 * (A + B) if lq > 0.0: t = 1.0 / (1.0 + exp(-lq)) else: e = exp(lq) t = e / (1.0 + e) return A + t * (B - A) cdef int _config3(double x, double fx, double xp, double fp, bint has2, double xp2, double fp2, double* out) except -1: cdef double xi, phi if not has2 or xp2 == xp or fp2 == fp: return 0 xi = (x - xp) / (xp2 - xp) phi = (fx - fp) / (fp2 - fp) if not (isfinite(xi) and isfinite(phi)): return 0 out[0] = xi out[1] = phi out[2] = fabs(x - xp) out[3] = fabs(fx - fp) return 1 cdef int _analogy_order(bint has_cur, double* cur, Buf* past, double* m_out) except -1: """1 and the order in m_out when the configuration is nearer to its own past than to the line phi = xi; the translation limit (m = inf) is on.""" cdef Py_ssize_t i, nb, ib = 0 cdef double xi, phi, d, dbest = 0.0, lx, lf cdef double* q nb = past.n // 4 if not has_cur or nb == 0: return 0 xi = cur[0] phi = cur[1] for i in range(nb): q = past.p + 4 * i d = _py_hypot(xi - q[0], phi - q[1]) if i == 0 or d < dbest: dbest = d ib = i q = past.p + 4 * ib if not ( _py_hypot(xi - q[0], phi - q[1]) < fabs(phi - xi) / sqrt(2.0)): return 0 if q[2] <= 0.0 or q[3] <= 0.0 or cur[2] <= 0.0 or cur[3] <= 0.0: return 0 lx = cur[2] / q[2] lf = cur[3] / q[3] if isfinite(lx) and isfinite(lf) and 0.0 < lf < 1.0 <= lx: m_out[0] = INFINITY return 1 if not (isfinite(lx) and isfinite(lf)) or lx == 1.0 or lf <= 0.0: return 0 m_out[0] = _plog(lf) / _plog(lx) return 1 cdef double _pattern_s(double step, double x, double xp, Buf* qs) except? -1.5: """The extrapolation share the history confirms (mode "pattern").""" cdef double d1 = x - xp, s_step, best = 0.0, mean cdef Py_ssize_t k, p, j cdef list sims if step == 0.0 or d1 == 0.0 or (step > 0.0) != (d1 > 0.0): return 0.0 s_step = _ratio_sim(fabs(step), fabs(d1)) k = qs.n for p in range(2, k // 2 + 1): sims = [_ratio_sim(qs.p[j], qs.p[j - p]) for j in range(p, k)] mean = _py_sum(sims) / len(sims) best = _pmax(best, mean) if best > 0.5: # max() can only keep or raise best, so the verdict is already fixed break return s_step if best > 0.5 else 0.0 cdef double _call(object f, double x, Py_ssize_t* n, Py_ssize_t limit) except? -1.5: if n[0] >= limit: raise _Budget() n[0] += 1 return f(x) # ---------------------------------------------------------------------------------------- # public interface # ---------------------------------------------------------------------------------------- def evaluation_bound(double a, double b, double xtol): """Worst-case number of evaluations of f guaranteed for the bracket [a, b] at tolerance xtol: twice bisection plus three, 2*ceil(log2(|b - a|/xtol)) + 3.""" cdef double w = fabs(b - a) if not (xtol > 0.0): raise ValueError("xtol must be positive") if w <= xtol: return 2 return 2 * int(ceil(log2(w / xtol))) + 3 def find_root(f, a, b, xtol=1e-12, maxfev=200): """Root of f in the bracket [a, b], f(a) and f(b) of opposite signs. Returns (x, nfev, converged): x the midpoint of the final bracket (or the exact root when f(x) == 0 was hit), nfev the number of evaluations of f, converged False only when the evaluation budget maxfev ran out first. Stops when the bracket is no wider than xtol; never stops on |f|. At most 2*ceil(log2(|b - a|/xtol)) + 3 evaluations.""" cdef double da = a, db = b, dxtol = xtol cdef Py_ssize_t limit = maxfev, n = 0 cdef double x = 0.5 * (da + db), fx, xp, fp, xp2 = 0.0, fp2 = 0.0 cdef bint has2 = False cdef double fa, fb, A, FA, B, FB, width, eps, allow, ratio, lh cdef long long half, n_half cdef double far, slope, sc, z_now, step = 0.0, mid, loc, brk, f_far, s_eff, longer, reach cdef double x_new, model_target = 0.0, mid_now, lin_mid, fnew, t_at, lo, hi, left, right cdef bint has_model, have_step, cur_ok, found, flag, steep, shallow, flag_over, enter cdef bint last_model = False, in_regime = False, analogy_seen = False cdef bint shallow_seen = False, has_target = False cdef double analogy_m = INFINITY, shallow_m = INFINITY, regime_target = 0.0, m_pair = 0.0 cdef double cur[4] cdef Buf configs, moves, qs configs.p = NULL; configs.n = 0; configs.cap = 0 moves.p = NULL; moves.n = 0; moves.cap = 0 qs.p = NULL; qs.n = 0; qs.cap = 0 try: try: fa = _call(f, da, &n, limit) fb = _call(f, db, &n, limit) if fa == 0.0: return da, n, True if fb == 0.0: return db, n, True if _same_sign(fa, fb): raise ValueError("f(a) and f(b) must straddle zero") A, FA, B, FB = da, fa, db, fb if fabs(fa) <= fabs(fb): x, fx, xp, fp = da, fa, db, fb else: x, fx, xp, fp = db, fb, da, fa # the projection guard: the bisection budget doubled, fixed by the initial bracket width = fabs(B - A) if dxtol == 0.0 or width <= 0.0: raise ValueError("the projection guard needs xtol and a positive initial width") eps = 0.5 * dxtol if width > dxtol: ratio = width / dxtol if not (ratio > 0.0) and not isnan(ratio): raise ValueError("math domain error") lh = ceil(log2(ratio)) if isnan(lh): raise ValueError("cannot convert float NaN to integer") if isinf(lh): raise OverflowError("cannot convert float infinity to integer") half = lh else: half = 0 n_half = ceil(2.0 * half) if n_half < half: n_half = half allow = pow(2.0, n_half) * eps while fabs(B - A) > dxtol: has_model = False far = B if fabs(B - x) >= fabs(A - x) else A slope = _slope3(x, fx, xp, fp, has2, xp2, fp2) sc = fabs(slope) * fabs(far - x) if sc == 0.0 or not isfinite(sc): x_new = 0.5 * (A + B) else: z_now = fabs(fx) / sc have_step = False mid = 0.5 * (A + B) # the flat regime, recognised by analogy with the run's own past cur_ok = _config3(x, fx, xp, fp, has2, xp2, fp2, cur) flag = (not in_regime) and last_model and fp != 0.0 and fx / fp >= 0.25 found = _analogy_order(cur_ok, cur, &configs, &m_pair) steep = found and m_pair >= 2.0 shallow = found and 0.0 < m_pair <= 2.0 / 3.0 flag_over = ((not in_regime) and last_model and fp != 0.0 and fx / fp <= -0.25) if steep: analogy_seen = True analogy_m = m_pair if shallow: shallow_seen = True shallow_m = m_pair enter = flag and analogy_seen if not enter and flag_over and shallow_seen: enter = True analogy_m = shallow_m if enter: in_regime = True if cur_ok: _push(&configs, cur[0]) _push(&configs, cur[1]) _push(&configs, cur[2]) _push(&configs, cur[3]) if in_regime: regime_target = _order_split(A, FA, B, FB, analogy_m) has_target = True step = regime_target - x have_step = True # the local picture contradicts the bracket if not have_step and x != xp: loc = (fx - fp) / (x - xp) brk = (FB - FA) / (B - A) if loc != 0.0 and (loc > 0.0) != (brk > 0.0): step = mid - x have_step = True # the half promised by the step against the half the bracket holds if not have_step: f_far = FB if far == B else FA if _sim(z_now) >= 0.5 and fabs(f_far) < fabs(fx): step = mid - x have_step = True if not have_step: step = (1.0 - _sim(z_now)) * (far - x) _push(&moves, x - xp) if moves.n >= 2 and moves.p[moves.n - 2] != 0.0 \ and moves.p[moves.n - 1] != 0.0: _push(&qs, fabs(moves.p[moves.n - 1]) / fabs(moves.p[moves.n - 2])) s_eff = _pattern_s(step, x, xp, &qs) if 0.0 < s_eff < 1.0: longer = step / (1.0 - s_eff) reach = fabs(far - x) step = longer if fabs(longer) <= reach else copysign(reach, step) if (fabs(step) < 0.5 * dxtol and last_model and fp != 0.0 and fabs(fx / fp) < 0.25): step = copysign(0.5 * dxtol, far - x) model_target = x + step has_model = True x_new = x + step # projection onto the interval the bisection budget still allows lo = _pmin(A, B) hi = _pmax(A, B) if isfinite(allow): left = hi - allow right = lo + allow if left > right: x_new = 0.5 * (lo + hi) else: x_new = _pmin(_pmax(x_new, left), right) if not isfinite(x_new) or not (_pmin(A, B) < x_new < _pmax(A, B)): x_new = 0.5 * (A + B) if x_new == x or x_new == xp: x_new = 0.5 * (A + B) if x_new == x or x_new == xp: return x_new, n, fabs(B - A) <= dxtol last_model = has_model and x_new == model_target mid_now = 0.5 * (A + B) lin_mid = 0.5 * (FA + FB) fnew = _call(f, x_new, &n, limit) if in_regime and has_target: t_at = (x_new - A) / (B - A) lin_mid = FA * (1.0 - t_at) + FB * t_at mid_now = x_new has_target = False if in_regime and x_new == mid_now and fabs(fnew - lin_mid) < fabs(fnew): # the point looks like a simple root: leave the regime in_regime = False analogy_seen = False allow *= 0.5 xp2, fp2, has2 = xp, fp, True xp, fp = x, fx x, fx = x_new, fnew if fnew == 0.0: return x_new, n, True if _opposite(FA, fnew): B, FB = x_new, fnew else: A, FA = x_new, fnew return 0.5 * (A + B), n, True except _Budget: return x, n, False finally: PyMem_Free(configs.p) PyMem_Free(moves.p) PyMem_Free(qs.p)