1.4.0 — the deprecation release: epistemic parity, zero breaks

EPIC3 Phase 2.  Every old name still works (removals are 2.0-only,
MIGRATION.md ships now):
- error-bar parity: density/cumulants/extreme grow with_err=True
  (per-probe scatter, atoms exact); kappa_w grows full=True (the whole
  per-direction curvature distribution)
- the trace split: trace_certified(f, support=) — one name per return shape
- defect.defect_barycentres (canonical; spectroscopy stays as legacy alias)
- resona.synthesize — the discoverable verb for from_measure
- lift.r_inverse / s_inverse — the R/S duals for spectral design
  (monotone-window check + bisection, round-trip < 1e-10)

108 tests; full-gallery ratchet vs baseline: zero metric diffs (timing
columns only) — all 1.4 additions are strictly additive.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: dimaq12

CHANGELOG.md

@@ -3,6 +3,39 @@
 All notable changes to resona.  The discipline throughout: every number below
 is printed by a test or a gallery stand, not asserted by hand.
 
+## [1.4.0] — 2026-06-12
+EPIC3 "the deprecation release" — honesty completion, zero breaks (every old
+name still works; removals happen in 2.0 only).  See MIGRATION.md.
+- **EPISTEMIC PARITY** — the 2.0 north star lands early: every stochastic
+  read now offers its scatter.  `density(xs, with_err=True)`,
+  `cumulants(order, with_err=True)` (nonlinearity propagated by per-probe
+  recomputation, not linearization), `extreme(with_err=True)` (a
+  reproducibility bar, stated as such), `wkernel.kappa_w(..., full=True)`
+  (the whole per-direction distribution — curvature anisotropy).
+- **The trace split**: `s.trace_certified(f, support=)` — one name per
+  return shape; `trace(certified=True)` keeps working through 1.x.
+- **`defect.defect_barycentres`** — the honest canonical name
+  (`spectroscopy` promised a spectrum it does not deliver; it remains as a
+  legacy alias until 2.0).
+- **`resona.synthesize`** — the discoverable verb for `from_measure`
+  (operator construction; both names stay).
+- **`lift.r_inverse` / `lift.s_inverse`** — the missing duals of the R/S
+  transforms for spectral design: monotone-window check + bisection to
+  machine tightness (round-trip < 1e-10 in tests).
+- **The measured retrofits** (EPIC3 Phase 0/1 — every delta measured in a
+  scratch BEFORE the gallery; falsified predictions recorded in
+  .audit/EPIC3_phase0.md): Chern marker via a precomputed Jackson/Chebyshev
+  projector (C = +0.985 at L=20 in seconds; 7× faster at L=12);
+  Heisenberg thermodynamics by deflating −H (ground states become exact
+  atoms: max Z err 19.4% → 0.09% at L=8, 16.3% → 0.67% at L=12); GP log-det
+  deflate=64 (0.84% → 0.53%); Koopman VdP at 4× data (|Δf| = 3e-4, finer
+  than the FFT reference's own 1e-3 bin).
+- **Persona landing pages** (docs/personas/): five half-page doors — quants,
+  physicists, ML, dynamics, numerical analysts — each mapping the persona's
+  first five tasks to five verified calls.
+- 108 tests; full-gallery ratchet: only the four retrofitted stands moved,
+  all by improvement; default paths bit-identical.
+
 ## [1.3.0] — 2026-06-12
 EPIC2 "the defect listens back" — every port stress-tested in the source
 corpus FIRST (FA/revise_stress/STRESS_REPORT.md: 4 of 5 gold survived; the

MIGRATION.md

@@ -16,8 +16,9 @@ bit-frozen default paths, same numbers.
 
 ## New surface in 1.4 (no old equivalent)
 
-- `resona.free.r_inverse` / `resona.free.s_inverse` — the missing duals of the
-  R/S-transforms, for spectral design (monotone, bisectable).
+- `resona.lift.r_inverse` / `resona.lift.s_inverse` — the missing duals of the
+  R/S-transforms, for spectral design (monotone window required, bisected to
+  machine tightness; they live beside `r_transform`/`s_transform`).
 - **Error-bar parity**: `density(xs, with_err=True)`, `cumulants(n, with_err=True)`,
   `extreme(with_err=True)`, `wkernel.kappa_w(..., full=True)` — every
   stochastic read now offers its scatter, the same way `trace`/`moment`/

docs/README.md

@@ -41,15 +41,17 @@ matvec = lambda v: A @ v          # your operator, however you can apply it
 | collapse the variance on a spiked operator (Hutch++) | `resona.of(mv, N, deflate=K)` | [reading-spectra](reading-spectra.md) | — |
 | high-resolution density without reorthogonalization | `resona.of(mv, N, k=256, engine="kpm")` | [reading-spectra](reading-spectra.md) | — |
 | get a CERTIFIED bracket (the answer provably inside) | `resona.quadform(mv, "inv", v, certified=True, support=(a,None))` | [precision-and-defects](precision-and-defects.md) | [`certified_logdet.py`](../examples/certified_logdet.py) |
-| certify the k-truncation of a trace estimate | `s.trace("log", certified=True, support=(a,None))` | [precision-and-defects](precision-and-defects.md) | [`certified_logdet.py`](../examples/certified_logdet.py) |
+| certify the k-truncation of a trace estimate | `s.trace_certified("log", support=(a,None))` | [precision-and-defects](precision-and-defects.md) | [`certified_logdet.py`](../examples/certified_logdet.py) |
+| error bars on density / cumulants / extremes | `s.density(xs, with_err=True)` etc. — every stochastic read has one | [reading-spectra](reading-spectra.md) | — |
 | resolve INTERIOR eigenvalues (spectrum slicing) | `s.zoom(a, b)` → polish nodes | [reading-spectra](reading-spectra.md) | [`spectra_to_machine_precision.py`](../examples/spectra_to_machine_precision.py) |
 | check how close a disorder computation is to the critical edge | `resona.subordination.contraction(s, xs, σ²)` → compare to 1 | [composing-operators](composing-operators.md) | — |
 | read the Koopman generator out of a legacy solver's error | `resona.defect.generator_read(P_n, P_2n, t, n)` | [precision-and-defects](precision-and-defects.md) | [`defect_spectroscopy.py`](../examples/defect_spectroscopy.py) |
-| read a band spectrum from the same error (blind-zone-free) | `resona.defect.spectroscopy(power, bands, coords)` | [precision-and-defects](precision-and-defects.md) | [`defect_spectroscopy.py`](../examples/defect_spectroscopy.py) |
+| read per-band barycentres from the same error (blind-zone-free) | `resona.defect.defect_barycentres(power, bands, coords)` | [precision-and-defects](precision-and-defects.md) | [`defect_spectroscopy.py`](../examples/defect_spectroscopy.py) |
 | follow eigenvalues through CROSSINGS along a parameter path | `resona.wkernel.track(A0, Bs, path)` | [inverse-problems](inverse-problems.md) | — |
 | size a safe parameter step (trust region for frozen W) | `resona.wkernel.kappa_w(A0, Bs, k0)` | [inverse-problems](inverse-problems.md) | — |
 | disorder-average a DOS (no realizations) | `resona.subordination.averaged_dos(sA,σ,xs)` | [composing-operators](composing-operators.md) | [`anderson_localization.py`](../examples/anderson_localization.py) |
-| CONSTRUCT an operator with a prescribed spectrum | inverse-CDF levels → `resona.from_measure(levels, 1/N)` | [inverse-problems](inverse-problems.md) | [`spectral_phenomena/operator_synthesis.py`](../examples/spectral_phenomena/operator_synthesis.py) |
+| CONSTRUCT an operator with a prescribed spectrum | inverse-CDF levels → `resona.synthesize(levels, 1/N)` | [inverse-problems](inverse-problems.md) | [`spectral_phenomena/operator_synthesis.py`](../examples/spectral_phenomena/operator_synthesis.py) |
+| invert an R/S-transform for spectral design | `resona.lift.r_inverse(s, target)` | [composing-operators](composing-operators.md) | — |
 | probe a NON-HERMITIAN operator (Markov, Koopman, damping) | `resona.cloud(mv, N)` → `.radius() .abscissa() .nodes` | [reading-spectra](reading-spectra.md) | [`science/koopman_dynamics.py`](../examples/science/koopman_dynamics.py) |
 | read a topological invariant matrix-free | P = `apply(H, step, v)` chains → Chern marker | [solving-and-evolving](solving-and-evolving.md) | [`quantum/chern_from_noise.py`](../examples/quantum/chern_from_noise.py) |
 | turn a TIME SERIES into an operator | `mv, rmv, r = resona.lift.koopman(snapshots)` | [lifting-nonlinear](lifting-nonlinear.md) | [`science/koopman_dynamics.py`](../examples/science/koopman_dynamics.py) |

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "resona"
-version = "1.3.0"
+version = "1.4.0"
 description = "The FFT of operators: probe, compose, read — matrix-free spectral computation with a built-in cost law."
 readme = "README.md"
 requires-python = ">=3.9"

resona/__init__.py

@@ -19,14 +19,16 @@
                                    # for self-adjoint operators; cloud for general)
     resona.solve.rayleigh_polish   # PRECISION — spend effort only on the defect
 """
-from .spectral import Spectral, apply, quadform, local_spectrum, local_density, from_measure, from_eigenbasis
+from .spectral import (Spectral, apply, quadform, local_spectrum, local_density,
+                       from_measure, from_eigenbasis, synthesize)
 from . import wkernel, lift, beta, defect, free, subordination, cost, flow, solve, thermal
 from .cloud import cloud, Cloud
 
 #: convenience: ``resona.of(matvec, N)`` == ``Spectral.of(matvec, N)``
 of = Spectral.of
 
-__version__ = "1.3.0"
+__version__ = "1.4.0"
 __all__ = ["Spectral", "of", "apply", "quadform", "local_spectrum", "local_density",
-           "from_measure", "from_eigenbasis", "cloud", "Cloud", "wkernel", "lift", "beta",
-           "defect", "free", "subordination", "cost", "flow", "solve", "thermal"]
+           "from_measure", "from_eigenbasis", "synthesize", "cloud", "Cloud",
+           "wkernel", "lift", "beta", "defect", "free", "subordination", "cost",
+           "flow", "solve", "thermal"]

resona/defect.py

@@ -181,12 +181,14 @@ def generator_read(P_n, P_2n, t, n, solver="be"):
     return (4.0 * n / t ** 2) * D
 
 
-def spectroscopy(power, bands, coords=None):
+def defect_barycentres(power, bands, coords=None):
     """COMPRESS each band of a defect power distribution to ONE coordinate —
-    its energy barycentre (the BDS read).  Despite the name, this does not
-    return a spectrum: it returns one (location, amplitude) pair per band —
-    that compression is exactly what survives noise where ratio estimators
-    die.
+    its energy barycentre (the BDS read).  One (location, amplitude) pair per
+    band — that compression is exactly what survives noise where ratio
+    estimators die.
+
+    (1.4) This is the canonical name; `spectroscopy` is the legacy alias —
+    it promised a spectrum it does not deliver — and is removed in 2.0.
 
     `power`  : |D̂|² — the defect's energy in the caller's diagonalizing
                basis (resona does NOT choose your basis: you transform).
@@ -219,3 +221,7 @@ def spectroscopy(power, bands, coords=None):
         kbar.append(float((coords[b] * w).sum() / tot))
         signal.append(float(np.sqrt(tot)))
     return np.array(kbar), np.array(signal)
+
+
+#: legacy alias (≤1.3 name) — same function; removed in 2.0 (see MIGRATION.md)
+spectroscopy = defect_barycentres

resona/lift.py

@@ -81,6 +81,52 @@ def s_transform(s, w):
     return float(out[0]) if np.isscalar(w) else out
 
 
+def r_inverse(s, value, w_max=1.0, n=2001):
+    """w such that R(w) = value — the dual of `r_transform`, for spectral
+    DESIGN (choose the coordinate w that realizes a target R).  Scalar or
+    array `value`.
+
+    Samples R on (0, w_max], requires monotonicity there (raises otherwise —
+    the inverse is ill-posed; shrink `w_max`), then bisects to machine
+    tightness.  R is analytic with R'(0) = κ₂ > 0, so a monotone window
+    always exists near 0."""
+    ws = np.linspace(w_max / n, w_max, n)
+    rs = r_transform(s, ws)
+    return _monotone_invert(lambda w: r_transform(s, w), ws, rs, value,
+                            "R", w_max)
+
+
+def s_inverse(s, value, w_max=1.0, n=2001):
+    """w such that S(w) = value — the dual of `s_transform` (positive
+    spectra), same contract as `r_inverse`."""
+    ws = np.linspace(w_max / n, w_max, n)
+    rs = s_transform(s, ws)
+    return _monotone_invert(lambda w: s_transform(s, w), ws, rs, value,
+                            "S", w_max)
+
+
+def _monotone_invert(fn, ws, rs, value, name, w_max):
+    d = np.diff(rs)
+    if not (np.all(d > 0) or np.all(d < 0)):
+        raise ValueError(f"{name} is not monotone on (0, {w_max}]: the "
+                         "inverse is ill-posed there — shrink w_max")
+    sign = 1.0 if d[0] > 0 else -1.0
+    v = np.atleast_1d(np.asarray(value, float))
+    lo_r, hi_r = min(rs[0], rs[-1]), max(rs[0], rs[-1])
+    if v.min() < lo_r or v.max() > hi_r:
+        raise ValueError(f"value outside {name}'s range [{lo_r:.6g}, "
+                         f"{hi_r:.6g}] on (0, {w_max}]")
+    idx = np.searchsorted(sign * rs, sign * v).clip(1, len(ws) - 1)
+    lo, hi = ws[idx - 1].astype(float), ws[idx].astype(float)
+    for _ in range(60):                                     # machine-tight
+        mid = 0.5 * (lo + hi)
+        below = sign * np.atleast_1d(fn(mid)) < sign * v
+        lo = np.where(below, mid, lo)
+        hi = np.where(below, hi, mid)
+    out = 0.5 * (lo + hi)
+    return float(out[0]) if np.isscalar(value) else out
+
+
 def free_convolution(sA, sB, order=6):
     """Moments of  A ⊞ B  (free additive convolution) from the two spectra ALONE.
 

resona/spectral.py

@@ -831,6 +831,8 @@ def trace(self, f, with_err: bool = False, certified: bool = False,
         plus ``support=(a, b)`` with the needed endpoint) → (lo, hi), the
         Gauss–Radau bracket (Golub–Meurant) of the K-TRUNCATION error: the
         fully-converged SLQ value of these same probes provably lies inside.
+        (1.4+: prefer `trace_certified` — one name per return shape; this
+        spelling remains through 1.x and is removed in 2.0.)
         It does NOT certify the Monte-Carlo probe scatter — that is a separate
         error source, reported by ``with_err``; the two are stated apart by
         design.  For an unconditional certificate of a probe-free quantity,
@@ -875,19 +877,69 @@ def trace(self, f, with_err: bool = False, certified: bool = False,
         stderr = float(np.std(ests, ddof=1) / np.sqrt(p))
         return total, stderr
 
+    def trace_certified(self, f, support=None):
+        """(lo, hi) — the Gauss–Radau bracket of Tr f(A)'s k-truncation.
+
+        The 1.4+ canonical name for ``trace(f, certified=True, support=…)``
+        (one name per return shape: `trace` returns a number, this returns
+        the bracket).  f must be a family name ('log', 'inv', 'sqrt',
+        'exp'); see `trace` for the full epistemics."""
+        return self.trace(f, certified=True, support=support)
+
     def moment(self, p: int, with_err: bool = False):
         """Tr(A^p).  ``with_err=True`` → (value, stderr), see `trace`."""
         return self.trace(lambda x: x ** p, with_err=with_err)
 
-    def density(self, xs, eta: float = 0.1) -> np.ndarray:
-        """Density of states ρ(x), Lorentzian-broadened by `eta`."""
+    def _probe_blocks(self):
+        """(sizes, n_atoms) — the per-probe node-block structure, validated."""
+        if not self.probe_sizes or len(self.probe_sizes) < 2:
+            raise ValueError("error bars need the probe structure from "
+                             "Spectral.of with probes >= 2")
+        used = int(np.sum(self.probe_sizes))
+        return list(self.probe_sizes), len(self.nodes) - used
+
+    def density(self, xs, eta: float = 0.1, with_err: bool = False):
+        """Density of states ρ(x), Lorentzian-broadened by `eta`.
+
+        ``with_err=True`` → (rho, stderr): the per-x standard error from the
+        independent-probe scatter (same epistemics as `trace`; deflate atoms
+        are exact and contribute no scatter)."""
         xs = np.atleast_1d(np.asarray(xs, float))
-        return (self.weights[None, :] * (eta / np.pi)
-                / ((xs[:, None] - self.nodes[None, :]) ** 2 + eta ** 2)).sum(1)
+        ker = (self.weights[None, :] * (eta / np.pi)
+               / ((xs[:, None] - self.nodes[None, :]) ** 2 + eta ** 2))
+        rho = ker.sum(1)
+        if not with_err:
+            return rho
+        sizes, _ = self._probe_blocks()
+        p = len(sizes)
+        ests, i = [], 0
+        for sz in sizes:
+            ests.append(p * ker[:, i:i + sz].sum(1))
+            i += sz
+        atom = ker[:, i:].sum(1)                              # exact, zero scatter
+        ests = np.stack([e + atom for e in ests])
+        return rho, np.std(ests, axis=0, ddof=1) / np.sqrt(p)
+
+    def extreme(self, with_err: bool = False):
+        """Extreme eigenvalues (Lanczos resolves these first / most reliably).
 
-    def extreme(self) -> tuple[float, float]:
-        """Extreme eigenvalues (Lanczos resolves these first / most reliably)."""
-        return float(self.nodes.min()), float(self.nodes.max())
+        ``with_err=True`` → ((lo, hi), (lo_err, hi_err)): the standard error
+        of the per-probe extreme reads — a REPRODUCIBILITY bar (how much the
+        read moves probe to probe), not a certified enclosure."""
+        lo, hi = float(self.nodes.min()), float(self.nodes.max())
+        if not with_err:
+            return lo, hi
+        sizes, _ = self._probe_blocks()
+        p = len(sizes)
+        atoms = self.nodes[int(np.sum(sizes)):]
+        mins, maxs, i = [], [], 0
+        for sz in sizes:
+            blk = self.nodes[i:i + sz]
+            cand = np.concatenate([blk, atoms]) if len(atoms) else blk
+            mins.append(float(cand.min())); maxs.append(float(cand.max()))
+            i += sz
+        return (lo, hi), (float(np.std(mins, ddof=1) / np.sqrt(p)),
+                          float(np.std(maxs, ddof=1) / np.sqrt(p)))
 
     # ── TRANSFORMS (the lifted coordinates, read off the measure) ─────────────
     def cauchy(self, z):
@@ -911,12 +963,31 @@ def s(self, w):
     r_transform = r          # full-name aliases (same methods; pick what reads
     s_transform = s          # better at the call site)
 
-    def cumulants(self, order: int = 6):
+    def cumulants(self, order: int = 6, with_err: bool = False):
         """Free cumulants κ_1..κ_order — the canonical coordinates in which
-        composition is ADDITION (κ_n(A⊞B) = κ_n(A) + κ_n(B))."""
+        composition is ADDITION (κ_n(A⊞B) = κ_n(A) + κ_n(B)).
+
+        ``with_err=True`` → (kappa, stderr): per-cumulant standard error from
+        the independent-probe scatter (each probe's moments → its cumulants;
+        the nonlinearity is propagated by recomputation, not linearization)."""
         from .free import free_cumulants
         m = [self.moment(p) / self.N for p in range(1, order + 1)]
-        return free_cumulants(m)
+        kappa = free_cumulants(m)
+        if not with_err:
+            return kappa
+        sizes, _ = self._probe_blocks()
+        p = len(sizes)
+        per_probe = []
+        for r in range(p):
+            m_r = []
+            for q in range(1, order + 1):
+                vals = self.weights * (self.nodes ** q)
+                i = int(np.sum(sizes[:r]))
+                blk = float(np.sum(vals[i:i + sizes[r]])) * p
+                atom = float(np.sum(vals[int(np.sum(sizes)):]))
+                m_r.append(blk + atom)
+            per_probe.append(free_cumulants(m_r))
+        return kappa, np.std(np.stack(per_probe), axis=0, ddof=1) / np.sqrt(p)
 
     # ── FLOW / DISORDER (the resolvent fixed point, vectorized) ───────────────
     def flow(self, t, xs, eta: float = 1e-3, g0=None):
@@ -1235,6 +1306,13 @@ def from_measure(nodes, weights, k=None):
     return _lanczos(lambda x: lam * x, np.sqrt(np.clip(w, 0.0, None)), k or len(lam))
 
 
+def synthesize(nodes, weights, k=None):
+    """CONSTRUCT an operator with the prescribed spectral measure — the
+    discoverable verb (1.4+) for `from_measure` (the same function; both
+    names stay, `from_measure` remains the precise synonym)."""
+    return from_measure(nodes, weights, k=k)
+
+
 def from_eigenbasis(eigenvalues, eigenvectors):
     """Reconstruct a Jacobi (tridiagonal) operator's band from its FULL
     eigendecomposition A = V·diag(λ)·Vᵀ — EXACT (machine precision) for ANY