"""Phase 7 — bootstrap robustness of the detected knee.
Rebuilds the curve as "its own smooth shape plus noise", resamples that noise
(IID residual bootstrap) and reruns the *full inner search* on each synthetic
curve. A knee is only robust if it is redetected almost every time, its
location has a tight interval, and the replicate knees are unimodal.
Author
------
Warith Harchaoui, <warithmetics@deraison.ai>
"""
from __future__ import annotations
from dataclasses import replace
from typing import List
import numpy as np
from .config import RobustKneeConfig
from .search import run_search
from .smoothing import signal_and_noise
from .types import BootstrapEvidence, PreparedCurve, Reason
def _greedy_1d_clusters(values: np.ndarray, tol: float) -> List[int]:
"""Sizes of greedy 1-D clusters (sorted values, running-median centers)."""
if values.size == 0:
return []
vs = np.sort(values)
sizes = []
current = [vs[0]]
center = vs[0]
for v in vs[1:]:
if abs(v - center) <= tol:
current.append(v)
center = float(np.median(current))
else:
sizes.append(len(current))
current = [v]
center = v
sizes.append(len(current))
return sorted(sizes, reverse=True)
[docs]
def bootstrap_knee(
prepared: PreparedCurve,
knee_x_norm: float,
config: RobustKneeConfig,
) -> BootstrapEvidence:
"""Assess knee stability under an IID residual bootstrap.
Parameters
----------
prepared : PreparedCurve
The observed normalized curve.
knee_x_norm : float
The accepted knee location. The replicate spread is reported
centred on it, and the median shift is measured against the knee of
the smoothed signal the replicates are drawn from.
config : RobustKneeConfig
Bootstrap thresholds and ``random_seed``.
Returns
-------
BootstrapEvidence
Detection rate, 90% interval, uni/multimodality rates, median shift,
and a pass flag with a reason code on failure.
"""
y = prepared.y_scaled
k = float(knee_x_norm)
yhat, residuals = signal_and_noise(y, config)
# The replicates are draws from "yhat + noise", so the location their
# median should be compared against is the knee *of yhat*, not the knee of
# the observed curve. The two differ by whatever the signal smoother did
# to the corner, and a smoother always moves it: see
# :func:`~elbow_helper.clustering._fine_scale_location`. Charging that
# displacement to the bootstrap reads as instability when nothing is
# unstable -- on the suite's 80-point knee-at-0.30 curve the replicates
# agree perfectly with each other (detection 1.00, primary cluster 1.00,
# CI width 0.063) yet sat 0.044 from the observed knee, over the 0.030
# tolerance, purely because the window-5 signal smoother had slid the
# corner 1.6 samples right. Referencing yhat's own knee measures what the
# gate is for: whether resampling noise *moves* the estimate.
reference = run_search(replace(prepared, y_scaled=yhat), config, confirm=False)
k_ref = (
float(reference.knee_x_norm)
if reference.detected and reference.knee_x_norm is not None
else k
)
rng = np.random.default_rng(config.random_seed)
knees: List[float] = []
b = config.bootstrap_replicates
for _ in range(b):
resampled = rng.choice(residuals, size=residuals.size, replace=True)
y_star = yhat + resampled
prepared_star = replace(prepared, y_scaled=y_star)
res = run_search(prepared_star, config, confirm=False)
if res.detected and res.knee_x_norm is not None:
knees.append(res.knee_x_norm)
detection_rate = len(knees) / b if b else 0.0
if not knees:
return BootstrapEvidence(
passes=False,
detection_rate=detection_rate,
ci90=(k, k),
ci90_width=1.0,
primary_cluster_rate=0.0,
secondary_cluster_rate=0.0,
median_shift=1.0,
knees=[],
reason=Reason.BOOTSTRAP_UNSTABLE,
)
knees_arr = np.array(knees)
lo, hi = np.percentile(knees_arr, [5, 95])
center = float(np.median(knees_arr))
ci90_width = float(hi - lo)
median_shift = float(abs(center - k_ref))
# Report the replicate *spread*, recentred on the knee it qualifies.
#
# The replicates re-find the knee of ``yhat``, so their percentiles are
# centred wherever the signal smoother put the corner -- at ``k_ref``, not
# at the ``k`` this interval is published alongside. Handing the raw
# percentiles out as ``ClearKnee.ci90`` therefore produced intervals that
# exclude the very location they are an interval for: on the measurement
# suite 20 of 28 accepted knees came back like ``knee_x=6.5,
# ci90=(7.5, 7.5)``. An interval that does not contain its own point
# estimate is not a statement a caller can act on -- it reads as two
# contradictory answers -- and the plotting module draws it as a band
# around a marker sitting outside the band.
#
# What the bootstrap actually measures is how far resampling the noise
# *moves* the estimate, which is the shape of the replicate distribution
# about its own centre. Translating that shape onto ``k`` keeps every
# measured quantity intact -- the width is unchanged to the last bit, so
# ``max_ci90_width`` still gates exactly what it gated -- while making
# containment automatic: ``lo <= center <= hi`` holds for any sample, so
# the translated interval brackets ``k`` always, with no clamping.
#
# The residual ``k - k_ref`` is not swept under the rug: it is precisely
# what ``median_shift`` measures, against its own threshold, just above.
ci90 = (k + (float(lo) - center), k + (float(hi) - center))
# Every threshold below bounds a distance in normalized x, which is
# quantized at one sample spacing; see RobustKneeConfig.positional.
n = prepared.n
sizes = _greedy_1d_clusters(
knees_arr, config.positional(config.cluster_tolerance, n)
)
primary_rate = sizes[0] / len(knees) if sizes else 0.0
secondary_rate = sizes[1] / len(knees) if len(sizes) > 1 else 0.0
detect_ok = detection_rate >= config.min_bootstrap_detection_rate
width_ok = ci90_width <= config.positional(config.max_ci90_width, n, samples=3.0)
primary_ok = primary_rate >= config.min_primary_cluster_rate
secondary_ok = secondary_rate <= config.max_secondary_cluster_rate
shift_ok = median_shift <= config.positional(config.max_bootstrap_median_shift, n)
reason = None
if not secondary_ok or not primary_ok:
reason = Reason.BOOTSTRAP_MULTIMODAL
elif not (detect_ok and width_ok and shift_ok):
reason = Reason.BOOTSTRAP_UNSTABLE
return BootstrapEvidence(
passes=detect_ok and width_ok and primary_ok and secondary_ok and shift_ok,
detection_rate=detection_rate,
ci90=ci90,
ci90_width=ci90_width,
primary_cluster_rate=primary_rate,
secondary_cluster_rate=secondary_rate,
median_shift=median_shift,
knees=knees,
reason=reason,
)