sim-numerics/restructure_article.py

"""Revision 1, step 3: the paper restructured and shortened.

Main text: abstract (<= 250 words), introduction, sections 2-5 (the base rule), 6 (the flat regime,
with the full listing of A'), 7 (out-of-sample evidence and the controls), 8 (conclusion), code
and data availability.  Appendices: A (where the evaluations go, where the similarity is read -
old sections 6 and 7), B (wall-clock - old 5.2), C (negative results - old 10), D (the
pre-registrations and their commits).  Sections move VERBATIM; cross-references and table
numbers are remapped in one pass; revision history is removed; new text is added only where the
revision asked for it.  Exact anchors, each found exactly once, or nothing is written.

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

import re

PATH = "article/sim_numerics.md"

# ------------------------------------------------------------------ 1. cross-references
SEC = {"5.3": "5.2", "7": "A.2", "7.1": "A.3", "7.2": "A.4", "7.3": "A.5", "7.4": "A.6",
       "8": "6", "8.1": "6.1", "8.2": "6.2", "8.3": "6.3", "8.4": "6.4", "8.5": "6.5", "8.6": "6.6",
       "9": "7", "9.1": "7.1", "10": "C", "11": "8"}
TAB = {"5": "A1", "6": "A2", "8": "5", "9": "6", "9a": "7", "9b": "8", "9c": "9"}


def remap(text):
    def sec(m):
        a, b = m.group(2), m.group(3)
        na = SEC.get(a, a)
        if na == "C":
            assert b is None
            return "Appendix C"
        out = m.group(1) + na
        if b:
            out += "–" + SEC.get(b, b)
        return out
    text = re.sub(r"(§§?)(\d+(?:\.\d+)?)(?:–(\d+(?:\.\d+)?))?", sec, text)

    def secword(m):
        na = SEC.get(m.group(1), m.group(1))
        return "Appendix C" if na == "C" else ("Appendix " + na if na.startswith("A") else "Section " + na)
    text = re.sub(r"Section (\d+(?:\.\d+)?)(?!\d)(?!\.\d)", secword, text)
    return re.sub(r"Table (9a|9b|9c|9|8|6|5)(?![0-9a-z])", lambda m: "Table " + TAB[m.group(1)], text)


# ------------------------------------------------------------------ 2. edits (new numbering)
EDITS = []


def E(old, new):
    EDITS.append((old, new))


# revision history and pronouns (revision: "him" -> the method)
E(" An earlier version gave each competitor the better of the library's implementation and ours, problem by problem. That concession backfired in the competitor's favour: our Chandrupatla stopped on a bracket twice as wide as everybody else's — the factor this protocol corrects for ITP in item 1 and had not corrected there — and the per-problem minimum let the looser run win. With the tolerance corrected our implementation and the library agree on all 199 problems exactly.",
  " Our own implementations of Chandrupatla's and Brent's methods agree with the libraries; the tables use the libraries.")
E(" For Brent the library is slightly dearer than our implementation; we use the library anyway, because it is the implementation a reader can check.", "")
E("Chandrupatla's own method and Dekker's are among the baselines; an earlier version of these tables was drawn up before either was added.",
  "Chandrupatla's own method and Dekker's are among the baselines.")
E("produced a variant that diverged from the original on 42 of the 199 problems and was therefore discarded rather than quoted; the question stays open. And an earlier version of this measurement estimated a per-evaluation overhead by subtracting the measured cost of $f$, which returned *negative* overheads for four of six methods — an impossibility that condemned the estimator, and the reason nothing is subtracted here.",
  "produced a variant that diverged from the original on 42 of the 199 problems and is not quoted. Nothing is subtracted: estimating a per-evaluation overhead by subtracting the measured cost of $f$ returns negative overheads for four of six methods.")
E("Chandrupatla does not win those families with his model — his traced run bisects",
  "Chandrupatla's method does not win those families with its model — its traced run bisects")
E("against his own 297, 233, 233", "against its own 297, 233, 233")
E("Chandrupatla halves on up to all of his steps", "Chandrupatla's method halves on up to all of its steps")
# the constant, stated once and honestly
E("which §3.1 declares as the method's one free constant.", "which §3.1 declares as the base rule's one numerical constant.")
E("**The multiplier is the method's one free numerical constant.** It is not \"the factor already claimed in the bound\" — that argument is circular, the bound is twice bisection because two was chosen — and measurement settles it (`measure_budget_factor.py`, Linux):",
  "**The multiplier is the base rule's one numerical constant**, and it is not neutral (`measure_budget_factor.py`, Linux):")
E("one declared constant and the better canonical result, or no constant, a tighter bound and the better wide-bracket result.",
  "one declared numerical constant and the better canonical result, or none, a tighter bound and the better wide-bracket result.")
# citations found missing on revision
E("Each of those answers carries constants: Illinois halves the stale ordinate,",
  "Each of those answers carries constants: Illinois [7] halves the stale ordinate,")
E("| 15 | Ridders | no |", "| 15 | Ridders [5] | no |")
E("| 237 | Illinois | yes |", "| 237 | Illinois [7] | yes |")
E("where Illinois, at the price of one constant", "where Illinois [7], at the price of one constant")
E("Two remarks apply to both forms.\n",
  "Three remarks apply to both forms.\n\n- **Guaranteed bounds of this kind are classical.** Bus and Dekker [21] gave bracketing methods whose "
  "worst case is four to five times bisection; the scheduled projection guarantees twice bisection plus one "
  "step, the earned one bisection plus one step.")
# the decision layer, counted honestly (revision point 2), and the full listing (point 4)
LISTING = r"""

Counting numerical constants does not measure the method's freedom, however. Beyond these numbers A′ is a decision layer of about ten discrete rules — the flag, the analogy and its memory, the entry, the exit, the overshoot branch, the split, the step pattern, the half, the contradiction, the tolerance step — chosen on the two development sets among the many candidates of Appendix C. That selection is why the method was tested once on a corpus fixed before its code was frozen, and against simpler rules put in its place (§7). Algorithm 2 is the complete method; `simroot/pure.py` and the compiled `simbracket` implement exactly this and agree bit for bit.

**Algorithm 2 (A′).** Notation of Algorithm 1; $S(z)=\max(0,1-|z|)$.

```
state: bracket [A,B], FA, FB;  points x, xp, xp2 and their values fx, fp, fp2
       memory: configs (list of configurations), moves (list of displacements)
       flags:  last_model, in_regime, analogy_seen (order m_a), shallow_seen (order m_s)
initialise as Algorithm 1; allow <- 2^n_max * xtol/2
while |B - A| > xtol:
    far <- end of [A,B] further from x;   mid <- (A+B)/2;   model <- false
    f'  <- derivative at x of the parabola through (x, xp, xp2)
    sigma <- |f'| * |far - x|
    if sigma is 0 or not finite:  xn <- mid
    else:
        z <- |fx| / sigma
        # the flat regime (6.1-6.3)
        c   <- (xi, phi, |x - xp|, |fx - fp|) of the last three points      (if defined)
        flag      <- not in_regime and last_model and fx/fp >= 1/4         unkept promise
        flag_over <- not in_regime and last_model and fx/fp <= -1/4        overshoot
        m <- order (6) of the analogous pair: the past configuration nearest to c
             (Euclidean in (xi, phi)), if it is nearer than the line phi = xi;
             m = inf when that pair shows L_f < 1 <= L_x (translation); else none
        if m >= 2:          analogy_seen <- true;  m_a <- m
        if 0 < m <= 2/3:    shallow_seen <- true;  m_s <- m
        if flag and analogy_seen:              in_regime <- true
        elif flag_over and shallow_seen:       in_regime <- true;  m_a <- m_s
        append c to configs
        if in_regime:                                   xn <- split point (7) with m = m_a
        elif x != xp and sign((fx-fp)/(x-xp)) != sign((FB-FA)/(B-A)):  xn <- mid
        elif S(z) >= 1/2 and |f(far)| < |fx|:           xn <- mid
        else:
            step <- (1 - S(z)) * (far - x)                                 eq. (4)
            append x - xp to moves
            if the ratios of successive move lengths repeat (the lag 2..half the history
               of highest mean ratio similarity has mean > 1/2), and step and x - xp
               point the same way:
                s <- min(|step|,|x-xp|) / max(|step|,|x-xp|)
                step <- step / (1 - s), at most |far - x| long
            if |step| < xtol/2 and last_model and |fx/fp| < 1/4:
                step <- sign(far - x) * xtol/2                             tolerance step
            xn <- x + step;  model <- true
        xn <- clamp(xn, hi - allow, lo + allow)                            projection (3)
    if xn outside (A,B) or xn in {x, xp}:  xn <- mid;  if that too is in {x, xp}: stop
    last_model <- model and xn is the model's own point
    fn <- f(xn)
    if in_regime:                                                          exit
        t <- (xn - A)/(B - A)
        if |fn - (FA(1-t) + FB t)| < |fn|:  in_regime <- false;  analogy_seen <- false
    allow <- allow/2;   (xp2, xp, x) <- (xp, x, xn);   keep the sign change in [A,B]
return midpoint(A,B)
```"""
E("it is verified with no violation on every corpus of this paper and under every safeguard of the grid.",
  "it is verified with no violation on every corpus of this paper and under every safeguard of the grid." + LISTING)
# related work named on revision, found and read
E("| the split (7) | false position on $\\operatorname{sign}(f)|f|^{1/m}$ [16, 17] |",
  "| the split (7) | false position on $\\operatorname{sign}(f)|f|^{1/m}$ [16, 17]; inside a bracket, with the multiplicity given, the SFRFm method of [22] |\n"
  "| the dimensionless configuration $(\\xi,\\varphi)$ of §6.2 | the ratios of Chandrupatla's acceptance band [3]; every scaled regula falsi written in two dimensionless parameters [22] |")
E("The methods for multiple roots we found are single-point iterations without a bracket, with the multiplicity assumed, estimated from derivatives, or removed by a transformation ([17];",
  "The methods for multiple roots we found are single-point iterations without a bracket, with the multiplicity assumed, estimated from derivatives, or removed by a transformation ([17, 23];")
E("Chandrupatla takes the decision to halve by an empirical acceptance band [3].",
  "inside a bracket, [22] needs the multiplicity given and recognises a multiple root only by the number of evaluations spent, against a third of bisection's; Chandrupatla takes the decision to halve by an empirical acceptance band [3].")
# the title
E("# Root Finding by Similarity and Analogy: A Bracketing Method with One Declared Constant",
  "# Root Finding by Similarity and Analogy: A Bracketing Method with a Data-Selected Decision Layer")

ABSTRACT = """Bisection places its next iterate at the midpoint because it uses nothing about $f$. We replace that rule by a similarity rule — the step grows with how dissimilar $f(x)$ is to the target 0, measured by a similarity $S$ with $S(0)=1$ — and derive its contraction factor in closed form. The derivation fixes a base rule before any benchmark is consulted: with the tent similarity and a three-point divided-difference scale it is a Newton step of Muller's order 1.839, clipped to the bracket and projected onto a budget of twice bisection whose multiplier is its one numerical constant. On the 154 problems of Alefeld, Potra and Shi it is level with TOMS 748 and 1.5% behind Chandrupatla's method; on wide brackets it is clearly behind. The gap is a decision, when to give the model up, and we take it by analogy: at a flat root the dimensionless configuration of the last three points recurs under zoom, at a simple root it tends to a line, and the analogous pair gives the order that says whether and where to split the bracket. This decision layer, about ten discrete rules selected on the two development sets, gives the method A′. On 240 problems fixed before its code was frozen it needs 3656 evaluations, against 5717 for Chandrupatla, 5535 for ModAB and 4488 for A′ with Lagouanelle's multiplicity estimate in place of the analogy. The lead lies in multiple and fractional roots, which two of that corpus's six families contain, and is small on simple roots. With the multiplicity known in advance, the classical transformation $\\operatorname{sign}(f)|f|^{1/m}$ is 2.4 to 6.8 times cheaper on multiple roots."""

INTRO = """## 1. Introduction

A bracketing root finder must decide, at every step, where inside the bracket to evaluate next. Bisection answers "in the middle", which is the correct answer when nothing is known. Every improvement since has answered "where a local model of $f$ says the root is": a secant in false position, an inverse quadratic in Brent [2], a cubic in TOMS 748 [1], an interpolate–truncate–project cycle in ITP [4], an Anderson–Björck step behind a linearity check in ModAB [10, 20]. Each of those answers carries constants: Illinois [7] halves the stale ordinate, Brent accepts an interpolation subject to three empirical inequalities, ITP is governed by $\\kappa_1,\\kappa_2,n_0$.

This paper asks a different question. Suppose the only thing available is a *similarity function* — a map $S$ from an outcome to $[0,1]$, where 1 means "identical to what was wanted" and 0 means "not similar at all", in the sense used for parameter-free control in [6]. Can a root finder be built out of $S$ alone, and how does it compare with methods built out of polynomial models and tuned constants?

The paper makes four claims, each measured.

1. A contraction law for the similarity step (§2) fixes a *base rule* — tent similarity, three-point scale, a projection safeguard — before any benchmark is consulted (§3); its order is Muller's.
2. Projecting the step onto the bisection budget, instead of replacing it, keeps the order and a worst case of twice bisection plus one step; the budget's multiplier is the base rule's one numerical constant, and a constant-free variant costs a measured price (§3.1, §5).
3. The base rule's gap on wide brackets is a decision — when to give the model up — and a similarity of the run to its own past takes it (§6). The resulting method, A′, is tested once on a corpus fixed before its code was frozen, against simpler detectors put in its place, and against the classical transformation with the multiplicity known in advance (§7).
4. Every table uses one criterion, a sign change within $10^{-6}$ of the answer, because the benchmarks' published roots score six problems unsolved for every method (§4, §A.5).

The search that produced these results — readings of the similarity that did not pay, the decomposition of the wide-bracket gap, wall-clock time, and the rejected rules — is reported in Appendices A–C. The code, the corpora and the pre-registrations are released (Appendix D and the availability statement)."""

AVAIL = """## Code and data availability

Everything behind this paper is released with it: the reference implementation of the base rule and of A′ (`simroot/pure.py`); `simbracket`, a compiled library of A′ (Cython; wheels for Linux and Windows, Python ≥ 3.10) whose results equal the reference bit for bit on all 868 problems of the six corpora; the corpora as code (`corpus_*.py`); every sweep script with its results file; the pre-registrations with the commits of Appendix D; and `verify_article.py`, which recomputes or re-reads every number stated in this paper and is part of the test suite. All numbers come from one Linux machine (§4). Repository: ⟨repository URL⟩."""

APPENDIX_D_HEAD = """## Appendix D. Pre-registrations and the frozen code

Every rule of §6 and every test of §7 was written down, with its criteria and stop rule, before its code and its run. Table D1 gives, for each document, the commit that first added it and the commit that first added the results its run wrote, both read from git by `check_prereg_order.py`: 25 of the 30 were committed before their results; four were committed together with them, so git cannot order them; one was decided on its smoke run and has no results file. The code of A′ was committed (`8027d7f`) after corpus 4 was defined (`657df10`) and before its results (`6d2531e`); `simroot/pure.py` did not change in between.

**Table D1.** Pre-registrations, in the order they were written.

| pre-registration | commit | results file | commit | order |
|---|---|---|---|---|
"""

REFS = """21. J. C. P. Bus and T. J. Dekker, "Two efficient algorithms with guaranteed convergence for finding a zero of a function," *ACM Trans. Math. Softw.*, vol. 1, no. 4, pp. 330–345, 1975.
22. J. M. Fernández-Díaz and C. O. Menéndez-Pérez, "A common framework for modified Regula Falsi methods and new methods of this kind," *Math. Comput. Simul.*, vol. 205, pp. 678–696, 2023.
23. A. Cordero, B. Neta and J. R. Torregrosa, "Memorizing Schröder's method as an efficient strategy for estimating roots of unknown multiplicity," *Mathematics*, vol. 9, 2570, 2021."""


def split_sections(text):
    """[(heading line or '' for the preamble, body)] in document order."""
    parts, cur_head, cur = [], "", []
    for ln in text.split("\n"):
        if ln.startswith("## ") or ln.startswith("### "):
            parts.append((cur_head, "\n".join(cur)))
            cur_head, cur = ln, []
        else:
            cur.append(ln)
    parts.append((cur_head, "\n".join(cur)))
    return parts


def clean(body):
    lines = body.strip("\n").split("\n")
    while lines and lines[-1].strip() in ("", "---"):
        lines.pop()
    while lines and lines[0].strip() in ("", "---"):
        lines.pop(0)
    return "\n".join(lines)


def main():
    src = open(PATH, encoding="utf8").read()
    src = remap(src)
    bad = [o[:80] for o, _ in EDITS if src.count(o) != 1]
    if bad:
        for b in bad:
            print("ANCHOR NOT FOUND EXACTLY ONCE:", repr(b))
        raise SystemExit(f"{len(bad)} anchors failed; nothing written")
    for o, n in EDITS:
        src = src.replace(o, n)
    # Table 1: caption, table, then the text about it
    cap = "$c_2$ the coefficient of $z^2$ when $c=0$.\n\n"
    t0 = "| $S(z)$ | $c$ | predicted | measured |"
    t1 = "\n\nThe last row is not a slow rate"
    i, j, k = src.index(cap) + len(cap), src.index(t0), src.index(t1)
    assert i < j < k and src.count(cap) == src.count(t0) == src.count(t1) == 1
    src = src[:i] + src[j:k] + "\n\n" + src[i:j].rstrip() + src[k:]

    secs = split_sections(src)
    by = {}
    for h, b in secs:
        key = re.match(r"#+ ([^ ]+)", h).group(1) if h else "PRE"
        assert key not in by, key
        by[key] = (h, clean(b))

    def S(key, heading=None):
        h, b = by[key]
        return (heading or h) + "\n\n" + b

    pre = by["PRE"][1]
    out = [pre, "## Abstract\n\n" + ABSTRACT, "---", INTRO, "---",
           S("2."), S("2.1"), S("2.2"), "---", S("3."), S("3.1"), "---", S("4."), "---",
           S("5."), S("5.1"), S("5.3", "### 5.2 The law across the whole grid"), "---",
           S("8.", "## 6. The flat regime, recognised by analogy")]
    for n in range(1, 7):
        h = by[f"8.{n}"][0].replace(f"### 8.{n} ", f"### 6.{n} ")
        out.append(S(f"8.{n}", h))
    out += ["---", S("9.", "## 7. Out-of-sample evidence"),
            S("9.1", by["9.1"][0].replace("### 9.1 ", "### 7.1 ")), "---",
            S("11.", "## 8. Conclusion"), S("Future"), "---", AVAIL, "---"]
    out += ["## Appendix A. Where the evaluations go, and where the similarity is read",
            S("6.", "### A.1 Where the evaluations go"),
            S("7.", "### A.2 Where the similarity is read")]
    for n, a in ((1, 3), (2, 4), (3, 5), (4, 6)):
        out.append(S(f"7.{n}", by[f"7.{n}"][0].replace(f"### 7.{n} ", f"### A.{a} ")))
    out += ["---", S("5.2", "## Appendix B. Wall-clock time, and why it is not the cost model"), "---",
            S("10.", "## Appendix C. Negative results"), "---"]
    rows = [ln for ln in open("results_prereg_order.txt", encoding="utf8").read().split("\n")
            if ln.startswith("| `")]
    assert len(rows) == 30
    out += [APPENDIX_D_HEAD + "\n".join(rows), "---", S("References") + "\n" + REFS]
    new = "\n\n".join(x.strip("\n") for x in out) + "\n"
    open(PATH, "w", encoding="utf8").write(new)
    print(f"written: {len(new.split())} words")


if __name__ == "__main__":
    main()