"""SIM root finding - the similarity-induced step, in the formulation of "Parameter-Free Similarity Control" (Sitnikov). Mapping of the article's machinery onto root finding ---------------------------------------------------- article here ------- ---- input x the target value y* = 0 ("the expectation") prototype c_k an evaluated value f(x_k) (value space) prototype memory y_k the abscissa x_k that produced it similarity s(x, c_k) similarity of f(x_k) to zero softmax weights w_k eq. (2) prediction sum w_k y_k the root estimate, eq. (3) So a root finder IS the article's controller run backwards: it is asked "which x is associated with the outcome 0?" and answers by normalized similarity aggregation over its memory of evaluated points. Two levels of prototype are implemented: POINT prototypes x_hat = sum_k w_k x_k Direct transcription of eq. (3). Exactness note: on a linear f this recovers the root exactly only for the inverse-distance similarity; every other similarity function biases it. Measured, not assumed. PAIR prototypes x_hat = sum_ij w_ij r_ij , r_ij = secant root of (i,j) Each bracketing PAIR of memory points is a local rule ("if the outcome changes sign between x_i and x_j then the root is near r_ij"), exactly parallel to a fuzzy rule table entry. This is the level at which the method departs from anything classical: regula falsi is the special case of a single pair, and the stale-endpoint failure of regula falsi disappears because a pair containing a large-|f| stale endpoint is by construction dissimilar to the target and is down-weighted. Parameter freedom ----------------- Similarity needs a scale in value space. It is taken from the data in play - the largest |f| currently in memory - never from a chosen constant. There is no sigma, no temperature, no threshold, no tuned coefficient anywhere below. """ from __future__ import annotations import math from .baselines import Counter, _Budget EPS = 2.220446049250313e-16 # -------------------------------------------------------------------------- # similarity functions on the value axis, all taking already-normalised # residuals z = f / scale and returning a similarity in a comparable range # -------------------------------------------------------------------------- def sim_neg_sq(z): """Negative squared Euclidean distance to the target, then softmax. This is the article's default (eq. 6) and its benchmark choice.""" return -z * z def sim_neg_abs(z): """Negative distance to the target, then softmax.""" return -abs(z) def sim_inv(z): """Inverse distance. Returned as a log-similarity so the caller's softmax reproduces plain inverse-distance weights. This is the unique member of the family that is exact on a linear f (see module docstring).""" return -math.log(abs(z) + EPS) def sim_inv_sq(z): return -2.0 * math.log(abs(z) + EPS) def sim_cosine(z): """Cosine similarity between the 2-vectors (1, f_k) and (1, 0), i.e. the article's eq. (4) applied to the value axis lifted to 2-D.""" return 1.0 / math.sqrt(1.0 + z * z) SIMILARITIES = { "neg_sq": sim_neg_sq, "neg_abs": sim_neg_abs, "inv": sim_inv, "inv_sq": sim_inv_sq, "cosine": sim_cosine, } def _softmax(scores): m = max(scores) ex = [math.exp(s - m) for s in scores] tot = sum(ex) if tot == 0.0 or not math.isfinite(tot): n = len(ex) return [1.0 / n] * n return [e / tot for e in ex] # -------------------------------------------------------------------------- # SIM root finder # -------------------------------------------------------------------------- def sim_root(f, a, b, xtol, maxfev, similarity="neg_sq", prototypes="pair", memory=None, guard=False): """Similarity-induced bracketed root finder. similarity : key into SIMILARITIES prototypes : "point" - eq. (3) verbatim, prototypes are evaluated points "pair" - rule-level prototypes, similarity on the value axis "pair2d" - rule-level prototypes, similarity on the value axis AND the abscissa axis (support tightness) memory : None = unbounded memory of every evaluated point; an int k = keep only the k most recent points. Unbounded is the default; the truncation is an ablation knob, not a tuned parameter. guard : enable the parameter-free safeguard. A similarity step is taken only while the bracket keeps halving; the moment a whole halving has failed to happen, one bisection step is forced. The only constant involved is 1/2, which is the definition of bisection rather than a tuned value, and it bounds the worst case at twice bisection. Compare: Brent, Ridders, TOMS748 and ITP all need a safeguard here too, and ITP's carries three tuning constants (kappa1, kappa2, n0). """ sfun = SIMILARITIES[similarity] two_d = prototypes == "pair2d" 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 # memory of every evaluated point mem = [(a, fa), (b, fb)] guard_width = abs(b - a) # width at the last acknowledged halving pending = False # a similarity step has been tried since then while abs(b - a) > xtol: lo, hi = (a, b) if a < b else (b, a) width = hi - lo if guard: if width <= 0.5 * guard_width: guard_width = width # the halving happened; nothing owed pending = False elif pending: # a similarity step was already spent without halving the # bracket -> pay it back with one bisection step x = 0.5 * (a + b) fx = g(x) mem.append((x, fx)) pending = False if fx == 0.0: return x, g.n, True if fa * fx < 0.0: b, fb = x, fx else: a, fa = x, fx continue pending = True # the memory is the WHOLE evaluation history, including points that # have fallen outside the current bracket: a stale point is not # deleted, it is merely dissimilar to the target and therefore # carries little weight. (Restricting the memory to the bracket # leaves exactly the two endpoints and collapses the method onto # regula falsi - verified in stage 1.) live = mem if memory is None else mem[-memory:] scale = max(abs(fi) for (_, fi) in live) if scale == 0.0: return live[0][0], g.n, True if prototypes == "point": scores = [sfun(fi / scale) for (_, fi) in live] w = _softmax(scores) x_new = sum(wi * xi for wi, (xi, _) in zip(w, live)) else: # pair prototypes cand, scores = [], [] for i in range(len(live)): xi, fi = live[i] for j in range(i + 1, len(live)): xj, fj = live[j] if fj == fi: continue r = (xi * fj - xj * fi) / (fj - fi) if not math.isfinite(r) or not (lo < r < hi): # a rule whose conclusion falls outside the proven # bracket is inadmissible, not merely unlikely continue cand.append(r) # a pair is similar to the target when BOTH of its # endpoints produced outcomes similar to the target sc = sfun(fi / scale) + sfun(fj / scale) if two_d: # ...and when its support is tight. Without this # term a flat tail (|f| ~ 1e-16 at 39 units from # the root) looks maximally similar to the target # and the method stalls there - measured, not # assumed. Both axes are normalised by the data # in play (largest |f| in memory; current bracket # width), exactly as the article normalises its # input space so that each axis contributes # equally. sc += sfun(abs(xi - xj) / width) scores.append(sc) if not cand: x_new = 0.5 * (a + b) else: w = _softmax(scores) x_new = sum(wi * ri for wi, ri in zip(w, cand)) if not math.isfinite(x_new) or not (lo < x_new < hi): x_new = 0.5 * (a + b) # never propose a point already in memory (would waste the budget) if any(x_new == xi for (xi, _) in mem): x_new = 0.5 * (a + b) if any(x_new == xi for (xi, _) in mem): # the bracket can no longer be split in floating point. # Report success only if it is actually inside the # tolerance; otherwise this is stagnation and must be # reported as a failure, not counted as a win. return x_new, g.n, abs(b - a) <= xtol x = x_new fx = g(x) mem.append((x, fx)) if fx == 0.0: return x, g.n, True if fa * fx < 0.0: b, fb = x, fx else: a, fa = x, fx return 0.5 * (a + b), g.n, True except _Budget: return x, g.n, False