sim-numerics/diag_contraction_similarity.py

"""Is the flat regime DETECTABLE by the similarity of two successive contraction ratios?

Valeri, 2026-09-21: "similarity should solve exactly this problem".  The similarity of a
triple to a linear configuration, |phi - xi|, cannot: fun1 on a wide bracket and the
pure powers x^9, x^19 look equally non-linear, yet a quadratic model works on the first
and fails on the second.  The flat regime shows in the HISTORY instead: with
rho_n = |f_n / f_{n-1}|, a superlinear run has rho falling from step to step (successive
ratios dissimilar), a linearly converging one has rho nearly constant (similar).  So the
candidate similarity for this place is

    S = min(rho_n, rho_{n-1}) / max(rho_n, rho_{n-1})        (the SIM ratio form)

Declared before the run, with no threshold anywhere:
  * AUC of S, flat states (fun3-fun6) against each other group: the probability that a
    random flat state has the higher S.  0.5 = no separation, 1.0 = perfect.
  * the threshold-free rule "the contraction did not speed up", rho_n >= rho_{n-1}: how
    often it fires in each group.  The bar to beat is the three-point order estimate,
    m >= 2: 0.2% of canonical states, 5.8% of flat ones (diag_flat_detectable.py).
Nothing in the method changes; states come from the behaviour-neutral trace hook.

    python -u diag_contraction_similarity.py
"""
from __future__ import annotations

from scipy.optimize._tstutils import _APS_TESTS_DICTS, _CHANDRUPATLA_TESTS_DICTS

from simroot.pure import pure_root

XTOL, MAXFEV, ACC = 1e-12, 200, 1e-6
CONFIGS = {
    "proposed (projection)": dict(similarity="ratio", scale="derived3", guard="project"),
    "reading (3), as for m>=2": dict(similarity="ratio", scale="derived3", guard="safe",
                                     prestep="slope", direction="dissimilar"),
}
FLAT = ("fun3", "fun4", "fun5", "fun6")


def states(f, a, b, r, cfg):
    """(S, contraction_did_not_speed_up) for every traced state with three nonzero f."""
    st = {"trace": []}
    x, _n, conv = pure_root(f, a, b, XTOL, MAXFEV, stats=st, **cfg)
    if not (conv and abs(float(x) - r) <= ACC):
        return []
    out = []
    for s in st["trace"]:
        if s["xp2"] is None or 0.0 in (s["fx"], s["fp"], s["fp2"]):
            continue
        rho_n = abs(s["fx"] / s["fp"])
        rho_p = abs(s["fp"] / s["fp2"])
        hi, lo = max(rho_n, rho_p), min(rho_n, rho_p)
        if hi == 0.0 or hi != hi:
            continue
        out.append((lo / hi, rho_n >= rho_p))
    return out


def auc(pos, neg):
    """P(random pos > random neg) + 0.5 P(tie), by ranks (Mann-Whitney)."""
    if not pos or not neg:
        return float("nan")
    allv = sorted([(v, 1) for v in pos] + [(v, 0) for v in neg])
    rank_sum, i = 0.0, 0
    while i < len(allv):
        j = i
        while j < len(allv) and allv[j][0] == allv[i][0]:
            j += 1
        mean_rank = (i + 1 + j) / 2.0
        rank_sum += mean_rank * sum(1 for k in range(i, j) if allv[k][1] == 1)
        i = j
    n1, n0 = len(pos), len(neg)
    return (rank_sum - n1 * (n1 + 1) / 2.0) / (n1 * n0)


def main():
    out = {}
    for cname, cfg in CONFIGS.items():
        groups = {"APS": [], "CHA fun1": [], "CHA flat fun3-6": [], "CHA other": []}
        for d in _APS_TESTS_DICTS:
            a, b = (float(v) for v in d["bracket"])
            f = (lambda ff, aa: lambda t: ff(t, *aa))(d["f"], d.get("args", ()))
            groups["APS"] += states(f, a, b, float(d["root"]), cfg)
        for d in _CHANDRUPATLA_TESTS_DICTS:
            a, b = (float(v) for v in d["bracket"])
            fam = str(d["ID"]).split(".")[0]
            g = ("CHA flat fun3-6" if fam in FLAT else "CHA fun1" if fam == "fun1"
                 else "CHA other")
            groups[g] += states(d["f"], a, b, float(d["root"]), cfg)
        print(f"\n== {cname}")
        print(f"   {'group':<18} {'states':>7} {'rho_n >= rho_(n-1)':>20}")
        for g, v in groups.items():
            n = max(len(v), 1)
            print(f"   {g:<18} {len(v):7d} {100.0 * sum(x[1] for x in v) / n:19.1f}%")
        flat = [x[0] for x in groups["CHA flat fun3-6"]]
        print("   AUC of S, flat states against:")
        res = {}
        for g in ("APS", "CHA fun1", "CHA other"):
            res[g] = auc(flat, [x[0] for x in groups[g]])
            print(f"     {g:<14} {res[g]:.3f}")
        out[cname] = {"groups": {g: (len(v), sum(x[1] for x in v)) for g, v in groups.items()},
                      "auc": res}
    print("\nbar to beat (three-point order estimate, m >= 2): 0.2% canonical, 5.8% flat")
    return out


if __name__ == "__main__":
    main()