sim-numerics/stage0_principle.py

"""Stage 0 - does the similarity-proportional step actually converge faster than bisection?

The principle under test (Valeri's formulation):
    in bisection the next point is the MIDPOINT of the bracket;
    in the similarity version the next point is placed proportionally to how
    SIMILAR the current result f(x) is to the expectation (zero):
        similar   -> small step
        dissimilar-> large step

Formalisation used here. Bracket [a,b] with f(a)*f(b) < 0. Define the
similarity of an endpoint's value to the target 0 in a SCALE-FREE way, using
only the data inside the bracket (no parameter):

    S(a) = |f(b)| / (|f(a)| + |f(b)|)      (a is similar to the root when |f(a)| is small)
    S(b) = |f(a)| / (|f(a)| + |f(b)|)      S(a) + S(b) = 1

Next point = similarity-weighted aggregation of the prototypes {a, b}:

    x = S(a)*a + S(b)*b

which is identical to stepping away from `a` by (1 - S(a)) of the bracket.
Bisection is the degenerate case S(a) = S(b) = 1/2, i.e. "no similarity
information used at all".

This file measures function evaluations to reach |x - root| <= tol, for
    - bisection          (S == 1/2, the null model)
    - sim-step           (the principle, raw)
and prints the per-problem numbers. No aggregate statistic is used to decide
anything; every problem is shown.
"""
from __future__ import annotations

from simroot.baselines import _opposite, _same_sign  # 2026-09-21 fix

import math
from dataclasses import dataclass
from typing import Callable

TOL = 1e-12
MAXFEV = 200


@dataclass
class Problem:
    name: str
    f: Callable[[float], float]
    a: float
    b: float
    root: float


PROBLEMS = [
    Problem("x^2 - 2                 ", lambda x: x * x - 2.0, 0.0, 2.0, math.sqrt(2.0)),
    Problem("cos(x) - x              ", lambda x: math.cos(x) - x, 0.0, 1.0, 0.7390851332151607),
    Problem("exp(x) - 2              ", lambda x: math.exp(x) - 2.0, 0.0, 2.0, math.log(2.0)),
    Problem("x^3 - x - 2             ", lambda x: x ** 3 - x - 2.0, 1.0, 2.0, 1.5213797068045676),
    Problem("sin(x) - 0.5            ", lambda x: math.sin(x) - 0.5, 0.0, 1.5, math.asin(0.5)),
    # the classical stalling case for regula falsi: strong convexity, root near
    # one end -> one endpoint never moves
    Problem("x^10 - 1   (STALL case) ", lambda x: x ** 10 - 1.0, 0.0, 1.3, 1.0),
    Problem("exp(-x)*(x-1) (flat end)", lambda x: math.exp(-x) * (x - 1.0), -1.0, 40.0, 1.0),
    Problem("atan(x)                 ", lambda x: math.atan(x), -1.0, 3.0, 0.0),
]


def _count(f):
    """Wrap f with an evaluation counter."""
    state = {"n": 0}

    def g(x):
        state["n"] += 1
        return f(x)

    return g, state


def bisection(f, a, b, tol=TOL, maxfev=MAXFEV):
    g, st = _count(f)
    fa, fb = g(a), g(b)
    if fa == 0.0:
        return a, st["n"]
    if fb == 0.0:
        return b, st["n"]
    x = 0.5 * (a + b)
    while st["n"] < maxfev:
        x = 0.5 * (a + b)
        fx = g(x)
        if fx == 0.0 or (b - a) <= tol:
            return x, st["n"]
        if _opposite(fa, fx):
            b, fb = x, fx
        else:
            a, fa = x, fx
    return x, st["n"]


def sim_step(f, a, b, tol=TOL, maxfev=MAXFEV):
    """The raw principle: x = S(a)*a + S(b)*b with S scale-free inside the bracket."""
    g, st = _count(f)
    fa, fb = g(a), g(b)
    if fa == 0.0:
        return a, st["n"]
    if fb == 0.0:
        return b, st["n"]
    x = 0.5 * (a + b)
    while st["n"] < maxfev:
        da, db = abs(fa), abs(fb)
        tot = da + db
        if tot == 0.0:
            return 0.5 * (a + b), st["n"]
        s_a = db / tot            # similarity of endpoint a to the root
        s_b = da / tot
        x_new = s_a * a + s_b * b
        # guard: keep the iterate strictly inside the bracket (numerical only,
        # not a tuning parameter)
        lo, hi = (a, b) if a < b else (b, a)
        if not (lo < x_new < hi):
            x_new = 0.5 * (a + b)
        x = x_new
        fx = g(x)
        if fx == 0.0 or abs(b - a) <= tol:
            return x, st["n"]
        if _opposite(fa, fx):
            b, fb = x, fx
        else:
            a, fa = x, fx
    return x, st["n"]


def run():
    print(f"{'problem':26s} {'bisect fev':>11s} {'bisect err':>11s} "
          f"{'sim fev':>9s} {'sim err':>11s}   verdict")
    print("-" * 90)
    wins = ties = losses = 0
    for p in PROBLEMS:
        xb, nb = bisection(p.f, p.a, p.b)
        xs, ns = sim_step(p.f, p.a, p.b)
        eb = abs(xb - p.root)
        es = abs(xs - p.root)
        # a method only "wins" if it reached the root at all
        ok_b = eb <= 1e-8
        ok_s = es <= 1e-8
        if ok_s and (not ok_b or ns < nb):
            verdict = "SIM faster"
            wins += 1
        elif ok_s and ok_b and ns == nb:
            verdict = "tie"
            ties += 1
        else:
            verdict = "SIM slower/failed"
            losses += 1
        print(f"{p.name:26s} {nb:11d} {eb:11.2e} {ns:9d} {es:11.2e}   {verdict}")
    print("-" * 90)
    print(f"SIM faster: {wins}   tie: {ties}   SIM slower/failed: {losses}"
          f"   (of {len(PROBLEMS)}; maxfev={MAXFEV})")


if __name__ == "__main__":
    run()