"""
Wavelet band decomposition (CWT + synchrosqueezed variant), ported near-verbatim from
c:\\DataS\\InstrumentSimulator\\backend\\app\\wavelet.py (project "Macro Causal Lab").

No project-specific dependencies — pure numpy/ssqueezepy, safe to call from any service.
"""
import logging
import warnings

import numpy as np
from ssqueezepy import cwt, icwt
from ssqueezepy.utils.cwt_utils import center_frequency

logger = logging.getLogger(__name__)

MIN_SCALES_PER_BAND = 4

# Fixed period (days) lower-bounds for bands 0..5, independent of the analysis
# window's length. Without this, dividing [scales.min(), scales.max()] into
# num_levels *equal* log-spaced buckets means "band i" covers a totally
# different real-world period range depending on how much history is fed in
# (a longer window resolves longer scales, which shifts every intermediate
# boundary upward) - so the same band index is not comparable across window
# lengths. Anchoring to absolute day-periods keeps "band i" meaning the same
# oscillation whether the window is 3 months or 3 years; only the last band's
# upper bound stays open, since a longer window can genuinely resolve slower
# cycles a short one cannot.
CANONICAL_PERIOD_EDGES_DAYS = [0.3, 1.5, 4.0, 12.0, 35.0, 90.0]


def _period_to_log_scale_fn(wavelet, scales, n, num_calib_points=16):
    """Build an interpolation from log(period_days) to log(scale) by sampling
    center_frequency at a handful of scales (not every scale - it's not
    vectorized and can be slow over hundreds/thousands of scales).
    """
    log_scales_full = np.log(scales)
    calib_log_scales = np.linspace(log_scales_full.min(), log_scales_full.max(), min(num_calib_points, len(scales)))
    calib_periods = np.array([1 / center_frequency(wavelet, scale=float(np.exp(s)), N=n) for s in calib_log_scales])
    calib_log_periods = np.log(calib_periods)

    def period_to_log_scale(period_days):
        return float(np.interp(np.log(period_days), calib_log_periods, calib_log_scales))

    return period_to_log_scale


def band_decompose(
    values: list[float],
    dates: list[str],
    num_levels: int = 4,
    wavelet: str = "gmw",
) -> dict:
    x = np.asarray(values, dtype=float)
    n = len(x)
    if n < 32:
        raise ValueError("Serie trop courte pour une analyse ondelette (32 points minimum).")

    x_mean = float(x.mean())
    xc = x - x_mean

    Wx, scales = cwt(xc, wavelet=wavelet)
    log_scales = np.log(scales)

    period_to_log_scale = _period_to_log_scale_fn(wavelet, scales, n)
    period_edges = list(CANONICAL_PERIOD_EDGES_DAYS[:num_levels])
    while len(period_edges) < num_levels:
        period_edges.append(period_edges[-1] * 2.5)

    # scales[0] = smallest scale = highest frequency = shortest period (impulses/noise);
    # scales[-1] = largest scale = lowest frequency = longest period (slow trend).
    bands = []
    for i in range(num_levels):
        lo = period_to_log_scale(period_edges[i])
        hi = period_to_log_scale(period_edges[i + 1]) if i + 1 < len(period_edges) else log_scales.max()
        mask = (log_scales >= lo) & (log_scales <= hi if i == num_levels - 1 else log_scales < hi)
        scales_band = scales[mask]

        reconstruction_failed = False
        if len(scales_band) >= MIN_SCALES_PER_BAND:
            try:
                recon = icwt(Wx[mask, :], wavelet=wavelet, scales=scales_band, x_len=n)
            except Exception:
                # Root-caused: icwt's internal scaletype inference recursively splits
                # `scales` looking for a log/linear transition point (infer_scaletype ->
                # logscale_transition_idx), and on a masked SUBSET of the original CWT
                # scales that recursive sub-slice can come out too short, crashing with
                # "attempt to get argmax of an empty sequence" instead of a clean error.
                # Confirmed by reproducing it directly against ssqueezepy 0.6.6: retrying
                # with a regenerated, exactly log-uniform array over the same range/count
                # sidesteps the fragile inference entirely (trivially satisfies its "is
                # this log-spaced" check) and reliably recovers a real reconstruction —
                # verified on 5 independently reproduced failures, all recovered.
                try:
                    clean_scales = np.geomspace(scales_band.min(), scales_band.max(), len(scales_band))
                    recon = icwt(Wx[mask, :], wavelet=wavelet, scales=clean_scales, x_len=n)
                    logger.warning(
                        "band_decompose: icwt failed for band %d with original scales (n=%d, wavelet=%s) — "
                        "recovered using a regenerated log-uniform scale array", i, len(scales_band), wavelet,
                    )
                except Exception:
                    logger.warning(
                        "band_decompose: icwt failed for band %d (n=%d scales, wavelet=%s) even after "
                        "retry — falling back to a zero band", i, len(scales_band), wavelet, exc_info=True,
                    )
                    recon = np.zeros(n)
                    reconstruction_failed = True
            scale_lo, scale_hi = float(scales_band.min()), float(scales_band.max())
        else:
            # Too few scales in this bucket for a stable reconstruction; report a
            # flat zero band (its content ends up folded into the residual) so the
            # band list always stays num_levels long.
            logger.warning(
                "band_decompose: only %d scale(s) in band %d (< %d minimum) — falling back to a zero band",
                len(scales_band), i, MIN_SCALES_PER_BAND,
            )
            recon = np.zeros(n)
            reconstruction_failed = True
            scale_lo, scale_hi = float(np.exp(lo)), float(np.exp(hi))

        freq_at_scale_hi = center_frequency(wavelet, scale=scale_hi, N=n)
        freq_at_scale_lo = center_frequency(wavelet, scale=scale_lo, N=n)
        period_low_days = round(1 / freq_at_scale_lo, 1) if freq_at_scale_lo else None
        period_high_days = round(1 / freq_at_scale_hi, 1) if freq_at_scale_hi else None
        label = (
            f"{period_low_days}-{period_high_days}j"
            if period_low_days is not None and period_high_days is not None
            else f"bande {i + 1}"
        )
        bands.append(
            {
                "index": i,
                "label": label,
                "period_low_days": period_low_days,
                "period_high_days": period_high_days,
                "series": [round(float(v), 6) for v in recon],
                "reconstruction_failed": reconstruction_failed,
            }
        )

    reconstructed_total = x_mean + sum(np.array(band["series"]) for band in bands)
    residual = x - reconstructed_total

    return {
        "dates": dates,
        "original": [round(float(v), 6) for v in x],
        "mean": round(x_mean, 6),
        "bands": bands,
        "residual": [round(float(v), 6) for v in residual],
        "wavelet": wavelet,
    }


def windowed_band_decompose(
    values: list[float],
    dates: list[str],
    window_size: int,
    num_levels: int = 4,
    wavelet: str = "gmw",
) -> dict:
    """Decompose a long series by running band_decompose independently on
    consecutive slices of `window_size` points, then concatenating the
    results. Each slice is analyzed on its own terms (own local mean, own
    CWT scale range), so precision near the reconstructed curve stays high
    even over a long history; the trade-off is a visible discontinuity at
    each slice boundary, which is acceptable for this diagnostic view.
    """
    n = len(values)
    window_size = max(32, min(window_size, n)) if window_size else n
    if n <= window_size:
        result = band_decompose(values, dates, num_levels, wavelet)
        result["window_size"] = window_size
        result["chunks"] = 1
        return result

    all_dates: list[str] = []
    all_original: list[float] = []
    all_residual: list[float] = []
    band_series: list[list[float]] = [[] for _ in range(num_levels)]
    band_meta: list[dict | None] = [None] * num_levels
    band_failed: list[bool] = [False] * num_levels

    chunk_count = 0
    start = 0
    while start < n:
        end = min(start + window_size, n)
        if 0 < n - end < 32:
            # Don't leave a too-short trailing slice: fold the remainder in.
            end = n
        chunk = band_decompose(values[start:end], dates[start:end], num_levels, wavelet)
        chunk_count += 1
        all_dates.extend(chunk["dates"])
        all_original.extend(chunk["original"])
        all_residual.extend(chunk["residual"])
        for band in chunk["bands"]:
            band_series[band["index"]].extend(band["series"])
            band_meta[band["index"]] = band
            # OR across chunks: if any slice couldn't reconstruct this band, flag the
            # whole concatenated series as unreliable rather than losing that signal
            # when a later chunk happens to succeed and overwrites band_meta.
            band_failed[band["index"]] = band_failed[band["index"]] or band.get("reconstruction_failed", False)
        start = end

    bands = []
    for i in range(num_levels):
        meta = band_meta[i] or {"label": f"bande {i + 1}", "period_low_days": None, "period_high_days": None}
        bands.append(
            {
                "index": i,
                "label": meta["label"],
                "period_low_days": meta["period_low_days"],
                "period_high_days": meta["period_high_days"],
                "series": band_series[i],
                "reconstruction_failed": band_failed[i],
            }
        )

    return {
        "dates": all_dates,
        "original": all_original,
        "mean": round(float(np.mean(all_original)), 6) if all_original else 0.0,
        "bands": bands,
        "residual": all_residual,
        "wavelet": wavelet,
        "window_size": window_size,
        "chunks": chunk_count,
    }


def rolling_causal_bands(
    values: list[float],
    dates: list[str],
    start_idx: int,
    lookback: int,
    num_levels: int = 4,
    wavelet: str = "gmw",
    step: int = 1,
) -> dict:
    """Walk-forward band decomposition with no look-ahead: the band value
    reported for day `t` is computed from a CWT run only on the trailing
    `lookback` points ending at `t` (never anything after `t`). This is
    deliberately much slower than `band_decompose` (one CWT per day instead
    of one for the whole range) - that cost is the point: a single
    whole-range CWT lets every point "see" the entire future through the
    transform's global/symmetric support, which makes any backtest built on
    it meaningless (the peaks/troughs it finds were computed with hindsight).
    `step` > 1 re-runs the CWT every `step` days and holds the last computed
    value in between, trading fidelity for speed on long ranges.
    """
    n = len(values)
    out_dates: list[str] = []
    out_original: list[float] = []
    band_series: list[list[float]] = [[] for _ in range(num_levels)]
    band_labels = [f"bande {i + 1}" for i in range(num_levels)]
    band_failed = [False] * num_levels
    last_tip: list[float] | None = None
    recomputations = 0

    for t in range(start_idx, n):
        window_start = t - lookback + 1
        if window_start < 0:
            continue  # not enough trailing history yet for a full lookback window
        if last_tip is None or (t - start_idx) % step == 0:
            window_values = values[window_start:t + 1]
            window_dates = dates[window_start:t + 1]
            result = band_decompose(window_values, window_dates, num_levels, wavelet)
            last_tip = [band["series"][-1] for band in result["bands"]]
            band_labels = [band["label"] for band in result["bands"]]
            for i, band in enumerate(result["bands"]):
                band_failed[i] = band_failed[i] or band.get("reconstruction_failed", False)
            recomputations += 1
        out_dates.append(dates[t])
        out_original.append(values[t])
        for i in range(num_levels):
            band_series[i].append(last_tip[i])

    bands = [
        {
            "index": i,
            "label": band_labels[i],
            "period_low_days": None,
            "period_high_days": None,
            "series": band_series[i],
            "reconstruction_failed": band_failed[i],
        }
        for i in range(num_levels)
    ]

    return {
        "dates": out_dates,
        "original": out_original,
        "bands": bands,
        "wavelet": wavelet,
        "lookback": lookback,
        "step": step,
        "recomputations": recomputations,
    }


# ---------------------------------------------------------------------------
# Synchrosqueezed variant (opt-in, method="ssq"). Everything above this line
# is completely untouched by what follows: band_decompose, windowed_band_decompose
# and rolling_causal_bands are the exact same functions used when method="cwt"
# (the default), so any existing simulation re-run with its original config
# hits the same code path and produces byte-identical results as before.
#
# Rationale: plain cwt/icwt smears a single true frequency's energy across
# several neighboring scales (that's *why* band_decompose needs canonical period
# edges and a minimum-scales guard - the boundaries are fuzzy). Synchrosqueezing
# reassigns that smeared energy onto its estimated true instantaneous frequency
# before splitting into bands, which should reduce cross-band leakage. It also
# exposes two things plain cwt does not: per-band energy |Tx|^2 (how dominant
# that cycle is *right now*, independent of direction) and a dominant-frequency
# ridge (which cycle length is currently winning, tracked over time).
#
# Separate period edges from the cwt path's CANONICAL_PERIOD_EDGES_DAYS: with
# fs=1.0 (one sample = one day, matching our daily data), ssq_cwt has a hard
# Nyquist floor at exactly 2.0 days (max resolvable frequency = fs/2). The cwt
# path's edges start at 0.3 days - below that floor, so reusing them here would
# leave the fastest band permanently empty. This does mean "band 0" covers a
# different absolute range under ssq than under cwt; that's an inherent
# consequence of ssq's stricter frequency floor, not a rounding choice.
# ---------------------------------------------------------------------------

CANONICAL_PERIOD_EDGES_DAYS_SSQ = [2.0, 4.0, 10.0, 25.0, 60.0, 120.0]


def band_decompose_ssq(
    values: list[float],
    dates: list[str],
    num_levels: int = 4,
    wavelet: str = "gmw",
) -> dict:
    from ssqueezepy import issq_cwt, ssq_cwt
    from ssqueezepy.ridge_extraction import extract_ridges

    x = np.asarray(values, dtype=float)
    n = len(x)
    if n < 32:
        raise ValueError("Serie trop courte pour une analyse ondelette (32 points minimum).")

    x_mean = float(x.mean())
    xc = x - x_mean

    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        # fs=1.0: one sample = one day, so ssq_freqs comes out directly in
        # cycles/day and 1/ssq_freqs is a period in days with no extra
        # calibration step (unlike band_decompose's scale->period conversion,
        # which needs _period_to_log_scale_fn precisely because plain cwt's
        # scale axis has no absolute unit by itself).
        Tx, _Wx, ssq_freqs, _scales = ssq_cwt(xc, wavelet=wavelet, fs=1.0)

    ssq_freqs = np.asarray(ssq_freqs)
    periods = 1.0 / ssq_freqs

    period_edges = list(CANONICAL_PERIOD_EDGES_DAYS_SSQ[:num_levels])
    while len(period_edges) < num_levels:
        period_edges.append(period_edges[-1] * 2.5)

    bands = []
    for i in range(num_levels):
        lo_period = period_edges[i]
        hi_period = period_edges[i + 1] if i + 1 < len(period_edges) else float(periods.max())
        mask = (periods >= lo_period) & (periods <= hi_period if i == num_levels - 1 else periods < hi_period)

        reconstruction_failed = False
        if mask.any():
            try:
                Tx_band = np.zeros_like(Tx)
                Tx_band[mask, :] = Tx[mask, :]
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore", RuntimeWarning)
                    recon = np.real(issq_cwt(Tx_band, wavelet=wavelet))
            except Exception:
                # Defensive, mirroring band_decompose: never let one band's edge
                # case fail the whole request — but log it and flag it (see
                # band_decompose's identical fix for why this used to be silent).
                logger.warning(
                    "band_decompose_ssq: issq_cwt failed for band %d (wavelet=%s) — "
                    "falling back to a zero band", i, wavelet, exc_info=True,
                )
                recon = np.zeros(n)
                reconstruction_failed = True
            energy = np.sum(np.abs(Tx[mask, :]) ** 2, axis=0)
        else:
            logger.warning("band_decompose_ssq: no frequencies fall in band %d — falling back to a zero band", i)
            recon = np.zeros(n)
            energy = np.zeros(n)
            reconstruction_failed = True

        period_low_days = round(float(lo_period), 1)
        period_high_days = round(float(hi_period), 1)
        bands.append(
            {
                "index": i,
                "label": f"{period_low_days}-{period_high_days}j",
                "period_low_days": period_low_days,
                "period_high_days": period_high_days,
                "series": [round(float(v), 6) for v in recon],
                "energy": [round(float(v), 8) for v in energy],
                "reconstruction_failed": reconstruction_failed,
            }
        )

    reconstructed_total = x_mean + sum(np.array(band["series"]) for band in bands)
    residual = x - reconstructed_total

    try:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)
            ridge_idxs = extract_ridges(Tx, ssq_freqs, n_ridges=1, bw=4)
        ridge_periods = 1.0 / ssq_freqs[ridge_idxs.ravel()]
        ridge_period_days = [round(float(v), 3) for v in ridge_periods]
    except Exception:
        ridge_period_days = [None] * n

    return {
        "dates": dates,
        "original": [round(float(v), 6) for v in x],
        "mean": round(x_mean, 6),
        "bands": bands,
        "residual": [round(float(v), 6) for v in residual],
        "wavelet": wavelet,
        "ridge_period_days": ridge_period_days,
    }


def rolling_causal_bands_ssq(
    values: list[float],
    dates: list[str],
    start_idx: int,
    lookback: int,
    num_levels: int = 4,
    wavelet: str = "gmw",
    step: int = 1,
) -> dict:
    """Synchrosqueezed counterpart to rolling_causal_bands - identical
    walk-forward, no-look-ahead structure (same day-by-day trailing window,
    same step logic), just calling band_decompose_ssq instead of
    band_decompose. rolling_causal_bands itself is untouched.
    """
    n = len(values)
    out_dates: list[str] = []
    out_original: list[float] = []
    band_series: list[list[float]] = [[] for _ in range(num_levels)]
    energy_series: list[list[float]] = [[] for _ in range(num_levels)]
    ridge_series: list[float | None] = []
    band_labels = [f"bande {i + 1}" for i in range(num_levels)]
    band_failed = [False] * num_levels
    last_tip: list[float] | None = None
    last_energy_tip: list[float] | None = None
    last_ridge_tip: float | None = None
    recomputations = 0

    for t in range(start_idx, n):
        window_start = t - lookback + 1
        if window_start < 0:
            continue
        if last_tip is None or (t - start_idx) % step == 0:
            window_values = values[window_start:t + 1]
            window_dates = dates[window_start:t + 1]
            result = band_decompose_ssq(window_values, window_dates, num_levels, wavelet)
            last_tip = [band["series"][-1] for band in result["bands"]]
            last_energy_tip = [band["energy"][-1] for band in result["bands"]]
            last_ridge_tip = result["ridge_period_days"][-1]
            band_labels = [band["label"] for band in result["bands"]]
            for i, band in enumerate(result["bands"]):
                band_failed[i] = band_failed[i] or band.get("reconstruction_failed", False)
            recomputations += 1
        out_dates.append(dates[t])
        out_original.append(values[t])
        for i in range(num_levels):
            band_series[i].append(last_tip[i])
            energy_series[i].append(last_energy_tip[i])
        ridge_series.append(last_ridge_tip)

    bands = [
        {
            "index": i,
            "label": band_labels[i],
            "period_low_days": None,
            "period_high_days": None,
            "series": band_series[i],
            "energy": energy_series[i],
            "reconstruction_failed": band_failed[i],
        }
        for i in range(num_levels)
    ]

    return {
        "dates": out_dates,
        "original": out_original,
        "bands": bands,
        "wavelet": wavelet,
        "lookback": lookback,
        "step": step,
        "recomputations": recomputations,
        "method": "ssq",
        "ridge_period_days": ridge_series,
    }


def _turning_points(series: list[float], smooth_days: int) -> list[tuple[int, int]]:
    """Indices (and new direction, +1/-1) where the sign of the `smooth_days`-lag slope of
    `series` flips. Using the raw day-to-day slope would flag noise as a "reversal"; lagging
    over a few days smooths that out."""
    turns: list[tuple[int, int]] = []
    prev_sign = 0
    for t in range(smooth_days, len(series)):
        slope = series[t] - series[t - smooth_days]
        sign = 1 if slope > 0 else (-1 if slope < 0 else 0)
        if sign != 0:
            if prev_sign != 0 and sign != prev_sign:
                turns.append((t, sign))
            prev_sign = sign
    return turns


def _average_cycle_days(turn_points: list[tuple[int, int]]) -> float | None:
    """Empirical average gap (in days) between consecutive SAME-direction turning points —
    peak-to-peak or trough-to-trough, i.e. one full oscillation — measured from the actual
    reconstructed series rather than assumed from the band's theoretical period_low/high
    labels (which describe the wavelet scale bucket, not necessarily the real cycle length
    this specific window's price action produced). None if there aren't at least 2
    same-direction turns to measure a gap from."""
    peaks = [t for t, d in turn_points if d == -1]    # a downturn = just passed a peak
    troughs = [t for t, d in turn_points if d == 1]   # an upturn = just passed a trough
    gaps = [seq[i + 1] - seq[i] for seq in (peaks, troughs) for i in range(len(seq) - 1)]
    return sum(gaps) / len(gaps) if gaps else None


def wavelet_reliability(
    values: list[float],
    dates: list[str],
    start_idx: int,
    lookback: int,
    num_levels: int = 4,
    wavelet: str = "gmw",
    method: str = "cwt",
    step: int = 1,
    smooth_days: int = 3,
    tolerance_pct: float = 0.10,
    min_confirm_horizon: int = 3,
) -> dict:
    """
    For every turning point a CAUSAL (walk-forward, no look-ahead) decomposition flags —
    i.e. what a live signal would have shown on that date, using only data available up to
    then — checks whether redoing the decomposition later still shows a same-direction
    reversal near the original date. The fraction confirmed is a per-band confidence score:
    CWT/SSQ reconstructions are well known to be least reliable right at the edge of the
    available data — exactly where a live "the band just turned" reading gets made.

    A fixed confirm horizon (e.g. always +10 days) is meaningless across bands whose natural
    oscillation period can range from ~1 day to several weeks: 10 days is many cycles for a
    fast band (trivially "confirmed", inflating its score) but less than half a cycle for a
    slow one (never has time to actually turn back, deflating its score) — the exact pattern
    that shows up as fast bands scoring ~100% and the slow band scoring ~10% for no real
    reliability reason. Instead, per band: measure its own historical average peak-to-peak
    (or trough-to-trough) period from the causal series, and use `avg_cycle * (1 +
    tolerance_pct)` as a MAXIMUM horizon — a same-direction reversal confirms the original one
    if it shows up ANYWHERE in the forward window [original date, original date + that
    horizon], not just near one specific point in it.
    """
    rolling_decomposer = rolling_causal_bands_ssq if method == "ssq" else rolling_causal_bands
    batch_decompose = band_decompose_ssq if method == "ssq" else band_decompose

    causal = rolling_decomposer(
        values, dates, start_idx=start_idx, lookback=lookback,
        num_levels=num_levels, wavelet=wavelet, step=step,
    )
    causal_dates = causal["dates"]
    date_to_global = {d: i for i, d in enumerate(dates)}

    bands_out = []
    for band in causal["bands"]:
        turn_points = _turning_points(band["series"], smooth_days)
        avg_cycle_days = _average_cycle_days(turn_points)

        if avg_cycle_days is None:
            bands_out.append({
                "label": band["label"], "n_tested": 0, "n_confirmed": 0, "confidence_pct": None,
                "avg_cycle_days": None, "confirm_horizon": None, "turns": [],
            })
            continue

        confirm_horizon = max(min_confirm_horizon, round(avg_cycle_days * (1 + tolerance_pct)))

        tested = []
        for t, direction in turn_points:
            d_date = causal_dates[t]
            g_idx = date_to_global.get(d_date)
            if g_idx is None or g_idx + confirm_horizon >= len(values):
                continue  # not enough real future data yet to judge this one

            end_idx = g_idx + confirm_horizon
            window_start = end_idx - lookback + 1
            if window_start < 0:
                continue
            window_values = values[window_start:end_idx + 1]
            window_dates = dates[window_start:end_idx + 1]
            try:
                hindsight = batch_decompose(window_values, window_dates, num_levels, wavelet)
            except ValueError:
                continue
            h_series = hindsight["bands"][band["index"]]["series"]

            # Forward-only window [original date, original date + horizon] — the reversal
            # can be confirmed anywhere in it, not just near a specific point. window_values
            # was built to end exactly at end_idx = g_idx + confirm_horizon, so h_series's
            # last index already IS "original date + horizon"; no separate bound needed.
            d_pos = g_idx - window_start
            hindsight_turns = _turning_points(h_series, smooth_days)
            matches = [idx for idx, new_dir in hindsight_turns if idx >= d_pos and new_dir == direction]
            confirmed = bool(matches)
            # Actual delay = how many days after the original date the confirming reversal
            # showed up (earliest match) — compared against avg_cycle_days below to check
            # whether "one cycle" is actually a well-calibrated horizon, or systematically
            # too long/short relative to how quickly a real reversal actually confirms.
            actual_delay = None
            height_diff = None
            if confirmed:
                confirm_idx = min(matches)
                actual_delay = confirm_idx - d_pos
                # Direction agreeing doesn't mean the peak/trough had the same amplitude —
                # the live (causal) reading can be a real reversal but distorted in height by
                # edge effects. Compare the band's own value at the original causal turning
                # point against its value at the confirming point in the hindsight
                # reconstruction (same units as the band, e.g. price terms).
                height_diff = abs(band["series"][t] - h_series[confirm_idx])
            tested.append({
                "date": d_date, "confirmed": confirmed,
                "actual_delay_days": actual_delay, "height_diff": height_diff,
            })

        n_tested = len(tested)
        n_confirmed = sum(1 for t in tested if t["confirmed"])
        delays = [t["actual_delay_days"] for t in tested if t["actual_delay_days"] is not None]
        avg_actual_delay = round(sum(delays) / len(delays), 1) if delays else None
        calibration_gap = round(sum(abs(avg_cycle_days - d) for d in delays) / len(delays), 1) if delays else None
        heights = [t["height_diff"] for t in tested if t["height_diff"] is not None]
        avg_height_diff = round(sum(heights) / len(heights), 6) if heights else None
        bands_out.append({
            "label": band["label"],
            "n_tested": n_tested,
            "n_confirmed": n_confirmed,
            "confidence_pct": round(n_confirmed / n_tested * 100, 1) if n_tested else None,
            "avg_cycle_days": round(avg_cycle_days, 1),
            "confirm_horizon": confirm_horizon,
            "avg_actual_delay_days": avg_actual_delay,
            "calibration_gap_days": calibration_gap,
            "avg_height_diff": avg_height_diff,
            "turns": tested,
        })

    return {
        "bands": bands_out,
        "smooth_days": smooth_days,
        "tolerance_pct": tolerance_pct,
        "method": method,
    }