#!/usr/bin/env python3
"""Offline, empirical modal voicing of complete stereo piano renders. See VOICING.md."""

import argparse
import csv
from dataclasses import dataclass
import hashlib
import io
import json
import math
import os
from pathlib import Path
import re
import struct
import sys
import tempfile

import numpy as np


KNOTS = (28, 68, 112)
REPRESENTATIVE_VELOCITIES = (13, 28, 35, 40, 45, 48, 53, 60, 68, 76, 84, 92, 100, 112, 116, 124)
PARTIALS = 240
TIMES = np.array([0.06, 0.35, 0.65, 1.0, 1.4])
WINDOW_SECONDS = np.array([0.12, 0.30, 0.30, 0.30, 0.30])
EARLY_SECONDS = 0.15
GAIN_LIMITS = (-36.0, 24.0)
GAIN_STEP = 12.0
DECAY_LIMITS = (0.1, 4.0)
DECAY_STEP = (0.25, 2.0)
PROBE_DECAY_STEP = (0.5, 2.0)
CSV_FIELDS = ("key", "partial", "pp_db", "mf_db", "ff_db", "decay_scale")
OUTPUTS = ("calibration.csv", "calibration_data.rs", "metadata.json", "summary.json")


def read_wav(path):
    """Read little-endian RIFF PCM16/24/32 or float32, including extensible WAV.

    Channels remain separate, as float64 in [-1, 1] for integer PCM. Float
    samples are not clipped. Reject damaged/unsupported files instead of guessing.
    """
    raw = Path(path).read_bytes()
    if len(raw) < 12 or raw[:4] != b"RIFF" or raw[8:12] != b"WAVE":
        raise ValueError(f"{path}: expected little-endian RIFF/WAVE")
    end = struct.unpack_from("<I", raw, 4)[0] + 8
    if end > len(raw):
        raise ValueError(f"{path}: truncated RIFF")
    fmt, payload = None, None
    pos = 12
    while pos + 8 <= end:
        tag, size = struct.unpack_from("<4sI", raw, pos)
        pos += 8
        if pos + size > end:
            raise ValueError(f"{path}: truncated WAV chunk")
        if tag == b"fmt ":
            fmt = raw[pos:pos + size]
        elif tag == b"data":
            if payload is not None:
                raise ValueError(f"{path}: multiple data chunks unsupported")
            payload = raw[pos:pos + size]
        pos += size + (size & 1)
    if fmt is None or len(fmt) < 16 or payload is None:
        raise ValueError(f"{path}: missing fmt/data chunk")
    encoding, channels, rate, byte_rate, align, bits = struct.unpack_from("<HHIIHH", fmt)
    if encoding == 0xFFFE:
        if len(fmt) < 40 or struct.unpack_from("<H", fmt, 16)[0] < 22:
            raise ValueError(f"{path}: truncated extensible format")
        valid_bits = struct.unpack_from("<H", fmt, 18)[0]
        if valid_bits not in (0, bits) or fmt[28:40] != bytes.fromhex("00001000800000aa00389b71"):
            raise ValueError(f"{path}: unsupported extensible PCM layout")
        encoding = struct.unpack_from("<I", fmt, 24)[0]
    if channels not in (1, 2) or rate <= 0 or bits not in (16, 24, 32):
        raise ValueError(f"{path}: unsupported channels/rate/bit depth")
    if align != channels * (bits // 8) or byte_rate != rate * align or len(payload) % align:
        raise ValueError(f"{path}: invalid PCM block alignment")
    if encoding == 3 and bits == 32:
        samples = np.frombuffer(payload, dtype="<f4").astype(np.float64)
    elif encoding == 1 and bits == 24:
        octets = np.frombuffer(payload, dtype=np.uint8).reshape(-1, 3).astype(np.int32)
        integers = octets[:, 0] | (octets[:, 1] << 8) | (octets[:, 2] << 16)
        samples = ((integers ^ 0x800000) - 0x800000).astype(np.float64) / 8388608.0
    elif encoding == 1 and bits in (16, 32):
        samples = np.frombuffer(payload, dtype=f"<i{bits // 8}").astype(np.float64) / 2 ** (bits - 1)
    else:
        raise ValueError(f"{path}: unsupported WAV encoding {encoding}/{bits}")
    samples = samples.reshape(-1, channels)
    if not samples.size or not np.isfinite(samples).all():
        raise ValueError(f"{path}: empty or non-finite audio")
    return samples, rate, {"encoding": encoding, "bits": bits, "channels": channels}


def onset_index(samples, rate):
    """First 1 ms stereo RMS block above -40 dB relative to first-0.5s peak."""
    block = max(1, round(rate * 0.001))
    count = min(len(samples), round(rate * 0.5)) // block
    if not count:
        raise ValueError("silent onset: audio shorter than one RMS block")
    power = np.mean(samples[:count * block].reshape(count, block, -1) ** 2, axis=(1, 2))
    if np.max(power) <= 1e-20:
        raise ValueError("silent onset: no signal in first 0.5 seconds")
    return int(np.flatnonzero(power > np.max(power) * 1e-4)[0] * block)


def power_spectrum(samples, rate, nfft=None):
    """One-sided stereo mean-square POWER per bin, with Parseval/Hann scaling."""
    size = len(samples)
    if size < 4:
        raise ValueError("FFT window too short")
    if nfft is None:
        nfft = 2 * (1 << (size - 1).bit_length())
    if nfft < size:
        raise ValueError("FFT size cannot truncate window")
    window = np.hanning(size)
    spectrum = np.fft.rfft(samples * window[:, None], n=nfft, axis=0)
    power = np.mean(np.abs(spectrum) ** 2, axis=1) / (nfft * np.sum(window ** 2))
    power[1:-1 if nfft % 2 == 0 else None] *= 2.0
    return np.fft.rfftfreq(nfft, 1.0 / rate), power


def midi_frequency(key):
    return 440.0 * 2.0 ** ((key - 69) / 12.0)


def partial_frequencies(f0, stiffness, count=PARTIALS):
    n = np.arange(1, count + 1, dtype=float)
    return f0 * n * np.sqrt(1.0 + stiffness * n * n)


def weighted_median(values, weights):
    order = np.argsort(values, kind="stable")
    values, weights = np.asarray(values)[order], np.asarray(weights)[order]
    cumulative = np.cumsum(weights)
    half = np.sum(weights) * 0.5
    tolerance = np.sum(weights) * 1e-12
    index = int(np.searchsorted(cumulative, half - tolerance))
    # Resolve an exactly balanced pair symmetrically, including roundoff. Otherwise
    # a constant dB shift can flip which of two slope estimates is the median.
    if index + 1 < len(values) and abs(cumulative[index] - half) <= tolerance:
        return float((values[index] + values[index + 1]) / 2)
    return float(values[index])


def robust_average(values, weights):
    values, weights = np.asarray(values), np.asarray(weights)
    center = weighted_median(values, weights)
    scale = max(0.05, 1.4826 * weighted_median(np.abs(values - center), weights))
    adjusted = weights * np.minimum(1.0, 1.5 * scale / np.maximum(np.abs(values - center), 1e-12))
    return float(np.sum(values * adjusted) / np.sum(adjusted))


def spectral_peaks(freq, power, duration):
    """Prominent local lines only; noise does not become a candidate harmonic."""
    candidates = np.flatnonzero((power[1:-1] > power[:-2]) & (power[1:-1] >= power[2:])) + 1
    candidates = candidates[power[candidates] > np.max(power) * 1e-6]
    found, strengths = [], []
    for index in candidates:
        # Exclude the whole close-unison cluster when estimating the floor.
        # A median through its center mistakes neighboring strings for noise.
        radius = max(14.0 / duration, freq[index] * 0.024)
        exclusion = max(4.0 / duration, freq[index] * 0.010)
        lo, hi = np.searchsorted(freq, [freq[index] - radius, freq[index] + radius])
        background = power[lo:hi][abs(freq[lo:hi] - freq[index]) > exclusion]
        if not len(background):
            continue
        floor = max(float(np.median(background)), np.max(power) * 1e-12, 1e-30)
        prominence = 10.0 * np.log10(power[index] / floor)
        if prominence < 15.0:
            continue
        logs = np.log(np.maximum(power[index - 1:index + 2], 1e-30))
        denominator = logs[0] - 2.0 * logs[1] + logs[2]
        delta = np.clip(0.5 * (logs[0] - logs[2]) / min(denominator, -1e-15), -0.5, 0.5)
        found.append(float(freq[index] + delta * (freq[1] - freq[0])))
        relative_db = 10.0 * np.log10(power[index] / np.max(power))
        strengths.append(float(np.clip((prominence - 15.0) / 15.0, 0, 1)
                               * np.clip((relative_db + 60.0) / 20.0, 0, 1)))
    return np.asarray(found), np.asarray(strengths)


def early_partial_lines(samples, rate, predicted, gaps):
    """Match early local lines to a fixed, independently fitted long-FFT curve."""
    segment = samples[:round(EARLY_SECONDS * rate)]
    freq, power = power_spectrum(segment, rate)
    duration = len(segment) / rate
    maximum = max(float(np.max(power)), 1e-30)
    peaks = np.flatnonzero((power[1:-1] > power[:-2]) & (power[1:-1] >= power[2:])) + 1
    centers, confidence = predicted.copy(), np.zeros_like(predicted)
    used = set()
    for p, center in enumerate(predicted):
        # Disjoint search bands and explicit peak ownership prevent aliases.
        tolerance = min(gaps[p] * 0.23, max(2.0 / duration, center * 0.003))
        radius = min(gaps[p] * 0.23, 2.5 * tolerance)
        if center + radius >= rate * 0.49:
            continue
        candidates = peaks[abs(freq[peaks] - center) <= radius]
        flank = (abs(freq - center) < gaps[p] * 0.45) & (abs(freq - center) > radius * 1.15)
        if not len(candidates) or np.count_nonzero(flank) < 4:
            continue
        floor = max(float(np.median(power[flank])), maximum * 1e-12, 1e-30)
        best = None
        for index in candidates:
            if int(index) in used:
                continue
            prominence = 10 * math.log10(max(power[index] / floor, 1e-30))
            relative = 10 * math.log10(max(power[index] / maximum, 1e-30))
            if prominence <= 15 or relative <= -60:
                continue
            logs = np.log(np.maximum(power[index - 1:index + 2], 1e-30))
            denominator = logs[0] - 2 * logs[1] + logs[2]
            delta = np.clip(0.5 * (logs[0] - logs[2]) / min(denominator, -1e-15), -0.5, 0.5)
            observed = float(freq[index] + delta * (freq[1] - freq[0]))
            if abs(observed - center) > radius:
                continue
            proximity = math.exp(-0.5 * ((observed - center) / tolerance) ** 2)
            quality = float(np.clip((prominence - 15) / 15, 0, 1)
                            * np.clip((relative + 60) / 20, 0, 1) * proximity)
            if best is None or quality > best[0]:
                best = quality, observed, int(index)
        if best is not None:
            confidence[p], centers[p], index = best
            used.add(index)
    return centers, confidence


def identify_partials(samples, rate, key):
    """Independent constrained stiff-string fit, then measured line neighborhoods.

    The stiff-string curve identifies harmonic numbers, not production frequencies.
    Local measured centers (including an independent p1) are retained for analysis.
    """
    segment = samples[round(0.1 * rate):round(1.2 * rate)]
    freq, power = power_spectrum(segment, rate)
    duration = len(segment) / rate
    peaks, strength = spectral_peaks(freq, power, duration)
    nominal = midi_frequency(key)
    n = np.arange(1, PARTIALS + 1, dtype=float)

    def match(predicted):
        if not len(peaks):
            return np.zeros(predicted.shape, dtype=int), np.zeros_like(predicted)
        right = np.clip(np.searchsorted(peaks, predicted), 0, len(peaks) - 1)
        left = np.maximum(right - 1, 0)
        index = np.where(abs(peaks[left] - predicted) < abs(peaks[right] - predicted), left, right)
        tolerance = np.maximum(2.0 / duration, predicted * 0.003)
        distance = abs(peaks[index] - predicted) / tolerance
        confidence = strength[index] * np.exp(-0.5 * distance ** 2) * (distance < 2.5)
        confidence *= predicted < rate * 0.48
        return index, confidence

    # Deterministic coarse grid avoids silently locking high harmonics to B=0.
    b_grid = np.r_[0.0, np.geomspace(1e-7, 0.01, 65)]
    f_grid = nominal * 2.0 ** (np.linspace(-45, 45, 37) / 1200.0)
    predicted = f_grid[:, None, None] * n * np.sqrt(1.0 + b_grid[None, :, None] * n * n)
    _, confidence = match(predicted)
    score = np.sum(confidence / np.sqrt(n), axis=-1)
    fi, bi = np.unravel_index(np.argmax(score), score.shape)
    f0, stiffness = float(f_grid[fi]), float(b_grid[bi])
    if not len(peaks):
        f0, stiffness = nominal, 0.0
    for _ in range(6):
        predicted = partial_frequencies(f0, stiffness)
        indices, confidence = match(predicted)
        valid = confidence > 0.15
        if np.count_nonzero(valid) < 3:
            break
        nv, measured = n[valid], peaks[indices[valid]]
        xscale = np.max(nv ** 2)
        design = np.column_stack([np.ones(len(nv)), nv ** 2 / xscale])
        weights = confidence[valid] / np.sqrt(nv)
        target = (measured / nv) ** 2
        for _ in range(4):
            root = np.sqrt(weights)
            coeff = np.linalg.lstsq(design * root[:, None], target * root, rcond=None)[0]
            fitted = design @ coeff
            residual = abs(target - fitted) / np.maximum(target, 1.0)
            weights = confidence[valid] / np.sqrt(nv) * np.minimum(1.0, 0.003 / np.maximum(residual, 1e-12))
        f0 = float(np.clip(math.sqrt(max(coeff[0], 1.0)), nominal * 2 ** (-45 / 1200), nominal * 2 ** (45 / 1200)))
        stiffness = float(np.clip(coeff[1] / xscale / (f0 * f0), 0.0, 0.01))
    predicted = partial_frequencies(f0, stiffness)
    indices, confidence = match(predicted)
    centers = predicted.copy()
    if len(peaks):
        centers[confidence > 0.1] = peaks[indices[confidence > 0.1]]
    supported = confidence > 0.25
    count = int(np.count_nonzero(supported))
    residual_cents = (1200 * np.log2(centers[supported] / predicted[supported]))
    error = float(np.sqrt(np.mean(residual_cents ** 2))) if count else 0.0
    fit_confidence = float(np.clip((count - 2) / 6, 0, 1) * math.exp(-(error / 8) ** 2))
    # With too few lines B is unidentified; do not extrapolate dubious upper modes.
    confidence[2:] *= 0.35 + 0.65 * fit_confidence
    gaps = np.minimum(np.diff(np.r_[0.0, predicted]), np.diff(np.r_[predicted, partial_frequencies(f0, stiffness, PARTIALS + 1)[-1]]))
    boundary = bool(stiffness >= 0.0099 or abs(1200 * math.log2(f0 / nominal)) >= 44.9)
    long_confidence = confidence.copy()
    early_centers, early_confidence = early_partial_lines(samples, rate, predicted, gaps)
    # Early peaks supply line evidence, never f0/B evidence. An unidentified or
    # boundary-constrained curve cannot assign new high partials reliably.
    early_confidence *= fit_confidence if not boundary and fit_confidence >= 0.25 else 0.0
    use_early = (confidence <= 0.25) & (early_confidence > confidence)
    centers[use_early] = early_centers[use_early]
    confidence[use_early] = early_confidence[use_early]
    return centers, confidence, gaps, {
        "f0_hz": f0, "B": stiffness, "confidence": fit_confidence,
        "supported_lines": count, "residual_cents_rms": error,
        "constraint_boundary": boundary,
        "long_line_confidence": long_confidence.tolist(),
        "early_line_confidence": early_confidence.tolist(),
        "early_line_used": use_early.tolist(),
    }


def robust_decay(times, power, confidence):
    """Positive empirical dB/s loss and confidence; reject rises, floors and beating."""
    good = (times >= 0.2) & (confidence > 0.2) & (power > 0)
    t, y, w = times[good], 10 * np.log10(np.maximum(power[good], 1e-30)), confidence[good]
    if len(t) < 3 or t[-1] - t[0] < 0.7:
        return 0.0, 0.0
    slopes, weights = [], []
    for i in range(len(t)):
        for j in range(i + 1, len(t)):
            slopes.append((y[j] - y[i]) / (t[j] - t[i]))
            weights.append(w[i] * w[j] * (t[j] - t[i]))
    slope = weighted_median(slopes, weights)
    intercept = weighted_median(y - slope * t, w)
    residual = abs(y - (intercept + slope * t))
    scatter = weighted_median(abs(np.asarray(slopes) - slope), weights)
    if slope >= -0.4 or np.max(np.diff(y)) > 2.0 or np.max(residual) > 3.0:
        return 0.0, 0.0
    quality = float(np.mean(w) * math.exp(-(float(np.sqrt(np.mean(residual ** 2))) / 1.5) ** 2)
                    / (1 + (scatter / (0.3 * abs(slope) + 0.3)) ** 2)
                    * np.clip((-slope - 0.4) / 1.5, 0, 1))
    return float(-slope), quality


@dataclass
class Measurement:
    power: np.ndarray
    confidence: np.ndarray
    power_upper_bound: np.ndarray
    location_confidence: np.ndarray
    loss: np.ndarray
    loss_confidence: np.ndarray
    times: np.ndarray
    window_samples: np.ndarray
    rate: int
    report: dict


def analyze(samples, rate, key):
    onset = onset_index(samples, rate)
    samples = samples[onset:]
    sizes = np.rint(WINDOW_SECONDS * rate).astype(int)
    starts = np.maximum(0, np.rint(TIMES * rate - sizes / 2).astype(int))
    if len(samples) < np.max(starts + sizes):
        raise ValueError("audio ends before final 1.4-second measurement window")
    centers, line_confidence, gaps, fit = identify_partials(samples, rate, key)
    predicted = partial_frequencies(fit["f0_hz"], fit["B"])
    location_confidence = np.zeros(PARTIALS)
    long_confidence = np.asarray(fit["long_line_confidence"])
    supported = np.flatnonzero(long_confidence > 0.25)
    if fit["confidence"] >= 0.25 and not fit["constraint_boundary"] and len(supported):
        # Interpolate only inside the independently identified harmonic range.
        location_confidence[supported[0]:supported[-1] + 1] = fit["confidence"]
        # A missing fundamental is located by at least three reliable low modes;
        # extrapolating down from p2..p8 is insensitive to high-mode stiffness.
        if np.count_nonzero(long_confidence[1:8] > 0.25) >= 3:
            location_confidence[0] = fit["confidence"]
        # An observed early high line supports its own location only; do not
        # extrapolate missing-reference ceilings through the intervening tail.
        early_supported = np.asarray(fit["early_line_used"]) & (line_confidence > 0.25)
        location_confidence[early_supported] = fit["confidence"] * line_confidence[early_supported]
    powers = np.zeros((len(TIMES), PARTIALS))
    confidences = np.zeros_like(powers)
    upper_bounds = np.zeros_like(powers)
    tracks = np.zeros_like(powers)
    times = np.zeros(len(TIMES))
    for ti, (start, size) in enumerate(zip(starts, sizes)):
        times[ti] = (start + (size - 1) / 2) / rate
        freq, power = power_spectrum(samples[start:start + size], rate)
        total = max(float(np.sum(power)), 1e-30)
        for p in range(PARTIALS):
            for bound, center in ((True, predicted[p]), (False, centers[p])):
                half_width = min(gaps[p] * 0.23, max(2.5 / (size / rate), center * 0.008))
                if (not bound and line_confidence[p] <= 0.05) or center + half_width >= rate * 0.49:
                    continue
                band = (freq >= center - half_width) & (freq <= center + half_width)
                flank = (abs(freq - center) < gaps[p] * 0.45) & (abs(freq - center) > half_width * 1.15)
                if np.count_nonzero(flank) < 4 or not np.any(band):
                    continue
                floor = max(float(np.median(power[flank])), total * 1e-12, 1e-30)
                noise = floor * np.count_nonzero(band)
                band_total = float(np.sum(power[band]))
                if bound:
                    # Keep all measured energy, including the floor, plus one
                    # extra local-floor allowance. Never interpret this as a line.
                    upper_bounds[ti, p] = band_total + noise
                    continue
                signal = max(band_total - noise, 0.0)
                snr = 10 * math.log10(max(signal / noise, 1e-30))
                relative = 10 * math.log10(max(signal / total, 1e-30))
                quality = np.clip((snr - 10) / 15, 0, 1) * np.clip((relative + 60) / 20, 0, 1)
                powers[ti, p] = signal
                confidences[ti, p] = quality * line_confidence[p]
                tracks[ti, p] = freq[band][np.argmax(power[band])]
    losses, loss_confidence = np.zeros(PARTIALS), np.zeros(PARTIALS)
    for p in range(PARTIALS):
        losses[p], loss_confidence[p] = robust_decay(times, powers[:, p], confidences[:, p])
    rms_slice = samples[round(0.05 * rate):round(0.45 * rate)]
    report = {"onset_samples": onset, "onset_seconds": onset / rate,
              "window_centers_seconds": times.tolist(), "window_samples": sizes.tolist(),
              "window_start_samples": starts.tolist(), "window_durations_seconds": (sizes / rate).tolist(),
              "pitch_fit": fit, "partial_centers_hz": centers.tolist(),
              "predicted_centers_hz": predicted.tolist(),
              "power_upper_bound": upper_bounds.tolist(),
              "location_confidence": location_confidence.tolist(),
              "tracked_peaks_hz": tracks.tolist(), "power": powers.tolist(),
              "power_confidence": confidences.tolist(), "loss_db_per_s": losses.tolist(),
              "loss_confidence": loss_confidence.tolist(),
              "early_stereo_rms": float(np.sqrt(np.mean(rms_slice ** 2)))}
    return Measurement(powers, confidences, upper_bounds, location_confidence,
                       losses, loss_confidence, times, sizes, rate, report)


def velocity_weights(velocities):
    velocities = np.asarray(velocities, dtype=float).reshape(-1)
    result = np.zeros((len(velocities), 3))
    for i, velocity in enumerate(velocities):
        if velocity <= KNOTS[0]:
            result[i, 0] = 1.0
        elif velocity >= KNOTS[-1]:
            result[i, -1] = 1.0
        else:
            lo = 0 if velocity < KNOTS[1] else 1
            fraction = (velocity - KNOTS[lo]) / (KNOTS[lo + 1] - KNOTS[lo])
            result[i, lo:lo + 2] = (1 - fraction, fraction)
    return result


def fit_velocity_gains(velocities, corrections, confidence):
    design = velocity_weights(velocities)
    confidence = np.asarray(confidence)
    if not np.any(confidence > 0):
        return np.zeros(3)
    # Small zero prior stabilizes weak observations and unsupported endpoint knots.
    support = design.T @ confidence
    ridge = 0.02 + 0.5 * np.maximum(0, 1 - support)
    matrix = design.T @ (confidence[:, None] * design) + np.diag(ridge)
    # No cross-knot prior: a loud layer must not fill an absent soft-layer partial.
    return np.linalg.solve(matrix, design.T @ (confidence * corrections))


def window_decay_bias(loss, measurement, ti):
    """Finite-Hann-window dB bias relative to power at the exact window center."""
    size = measurement.window_samples[ti]
    t = (np.arange(size) - (size - 1) / 2) / measurement.rate
    w2 = np.hanning(size) ** 2
    exponent = -math.log(10) / 10 * loss * t
    shift = float(np.max(exponent))
    return 10 / math.log(10) * (shift + math.log(float(np.sum(w2 * np.exp(exponent - shift)) / np.sum(w2))))


def validate_measurement_windows(measurements):
    baseline = measurements[0]
    for measurement in measurements:
        times, sizes = np.asarray(measurement.times), np.asarray(measurement.window_samples)
        if (times.ndim != 1 or len(times) < 2 or sizes.shape != times.shape
                or not np.isfinite(times).all() or np.any(np.diff(times) <= 0)
                or not np.isfinite(sizes).all() or np.any(sizes < 4) or np.any(sizes != np.rint(sizes))
                or measurement.rate <= 0
                or any(np.shape(array) != (len(times), PARTIALS) for array in
                       (measurement.power, measurement.confidence, measurement.power_upper_bound))):
            raise ValueError("measurements have invalid time/window shapes")
        if not (np.array_equal(times, baseline.times)
                and np.array_equal(sizes / measurement.rate, baseline.window_samples / baseline.rate)):
            raise ValueError("measurements must use the same time windows and durations")


def fit_note(pairs, reference_offset_db, previous_gain=None, previous_decay=None):
    """Fit residual strike dB and one shared empirical decay ratio per partial."""
    for _, model, reference in pairs:
        validate_measurement_windows((model, reference))
    previous_gain = np.zeros((3, PARTIALS)) if previous_gain is None else previous_gain
    previous_decay = np.ones(PARTIALS) if previous_decay is None else previous_decay
    log_updates, decay_confidence = np.zeros(PARTIALS), np.zeros(PARTIALS)
    for p in range(PARTIALS):
        ratios, weights = [], []
        for _, model, reference in pairs:
            weight = model.loss_confidence[p] * reference.loss_confidence[p]
            if weight > 0.05 and model.loss[p] > 0 and reference.loss[p] > 0:
                ratios.append(math.log(reference.loss[p] / model.loss[p]))
                weights.append(weight)
        if weights:
            center = robust_average(ratios, weights)
            dispersion = weighted_median(abs(np.asarray(ratios) - center), weights)
            quality = min(1.0, sum(weights)) / (1 + (dispersion / 0.35) ** 2)
            log_updates[p] = center * quality
            decay_confidence[p] = quality
    # Light smoothing only on already supported, low-confidence decay estimates.
    regularized = log_updates.copy()
    for p in range(1, PARTIALS - 1):
        if 0 < decay_confidence[p] < 0.8:
            neighbors = decay_confidence[p - 1:p + 2:2]
            if np.sum(neighbors) > 0:
                amount = 0.15 * (1 - decay_confidence[p])
                regularized[p] = (1 - amount) * log_updates[p] + amount * np.average(log_updates[p - 1:p + 2:2], weights=neighbors)
    decay_step = np.clip(np.exp(regularized), *DECAY_STEP)
    decay = np.clip(previous_decay * decay_step, *DECAY_LIMITS)
    applied_ratio = decay / previous_decay
    gain = np.zeros((3, PARTIALS))
    partial_report = []
    for p in range(PARTIALS):
        velocities, corrections, weights = [], [], []
        cut_corrections, cut_weights = [], []
        for velocity, model, reference in pairs:
            values, quality = [], []
            cut_values, cut_quality = [], []
            ratio = applied_ratio[p]
            can_compensate = model.loss_confidence[p] > 0.05
            for ti in (0, 1):
                weight = model.confidence[ti, p] * reference.confidence[ti, p]
                if abs(ratio - 1) > 1e-6 and not can_compensate:
                    continue
                delta = 0.0
                if can_compensate:
                    delta = (-model.loss[p] * (ratio - 1) * model.times[ti]
                             + window_decay_bias(model.loss[p] * ratio, model, ti)
                             - window_decay_bias(model.loss[p], model, ti))
                if reference.confidence[ti, p] <= 0.2:
                    # Censored reference energy can only support attenuation of
                    # a measured model line at a reliably located reference band.
                    upper = reference.power_upper_bound[ti, p]
                    if model.confidence[ti, p] > 0.25 and reference.location_confidence[p] >= 0.25 and upper > 0:
                        ceiling = 10 * math.log10(upper / model.power[ti, p]) + reference_offset_db
                        if ceiling < 0 and ceiling - delta < 0:
                            cut_values.append(ceiling - delta)
                            cut_quality.append(0.25 * model.confidence[ti, p] * reference.location_confidence[p])
                elif weight > 0:
                    correction = 10 * math.log10(reference.power[ti, p] / model.power[ti, p]) + reference_offset_db - delta
                    values.append(correction)
                    quality.append(weight)
            velocities.append(velocity)
            if quality:
                corrections.append(robust_average(values, quality))
                weights.append(float(np.mean(quality)))
            else:
                corrections.append(0.0)
                weights.append(0.0)
            # Choose the least restrictive eligible early-window bound.
            # These are never decay evidence.
            cut_corrections.append(max(cut_values) if cut_values else 0.0)
            cut_weights.append(float(np.mean(cut_quality)) if cut_quality else 0.0)
        velocity_fit = fit_velocity_gains(velocities + velocities,
                                          np.asarray(corrections + cut_corrections), weights + cut_weights)
        # An intermediate layer can constrain two knots algebraically. Do not
        # let it boost spectral noise at a knot whose own render lacks the mode.
        knot_evidence = np.array([weights[velocities.index(v)] if v in velocities else 0.0 for v in KNOTS])
        boost_guard = np.clip(knot_evidence / 0.25, 0, 1)
        raw_update = np.where(velocity_fit > 0, velocity_fit * boost_guard, velocity_fit)
        # Every knot touching a censored observation must remain cut-only. This
        # also guards intermediate velocities against boosts from other layers.
        design = velocity_weights(velocities)
        cut_knots = design.T @ np.asarray(cut_weights) > 0
        raw_update[cut_knots] = np.minimum(raw_update[cut_knots], 0.0)
        update = np.clip(raw_update, -GAIN_STEP, GAIN_STEP)
        gain[:, p] = np.clip(previous_gain[:, p] + update, *GAIN_LIMITS)
        applied_by_velocity = design @ (gain[:, p] - previous_gain[:, p])
        partial_report.append({
            "partial": p + 1, "velocities": velocities,
            "residual_db_by_velocity": corrections, "gain_confidence_by_velocity": weights,
            "cut_only_upper_db_by_velocity": cut_corrections,
            "cut_only_confidence_by_velocity": cut_weights,
            "cut_only_residual_db_before": np.maximum(0, -np.asarray(cut_corrections)).tolist(),
            "cut_only_residual_db_after": np.where(np.asarray(cut_weights) > 0,
                np.maximum(0, applied_by_velocity - cut_corrections), 0).tolist(),
            "knot_cut_only": cut_knots.tolist(),
            "velocity_fit_update_db": velocity_fit.tolist(),
            "knot_boost_evidence_factors": boost_guard.tolist(),
            "unclamped_gain_update_db": raw_update.tolist(),
            "applied_gain_update_db": (gain[:, p] - previous_gain[:, p]).tolist(),
            "gain_clamped": bool(np.any(abs(raw_update - (gain[:, p] - previous_gain[:, p])) > 1e-8)),
            "decay_confidence": float(decay_confidence[p]),
            "unclamped_decay_ratio": float(np.exp(regularized[p])),
            "applied_decay_ratio": float(applied_ratio[p]),
            "decay_clamped": bool(abs(regularized[p] - math.log(applied_ratio[p])) > 1e-8),
        })
    return gain, decay, partial_report


def response_lstsq(design, target, weight, prior):
    """Small ridge fit with 3 dB Huber residual weights (no acoustic model)."""
    effective = weight.copy()
    update = np.zeros(design.shape[1])
    for _ in range(12):
        update = np.linalg.solve(design.T @ (effective[:, None] * design) + np.diag(prior),
                                 design.T @ (effective * target))
        effective = weight * np.minimum(1, 3 / np.maximum(abs(target - design @ update), 1e-12))
    return update, effective


def response_condition(derivative, weight, times):
    """Center within each velocity, so layer differences cannot identify decay."""
    centered = np.zeros_like(derivative)
    eligible = np.zeros_like(weight)
    spans = []
    for i, row in enumerate(weight):
        good = row > 0.05
        span = float(np.ptp(times[i, good])) if np.any(good) else 0.0
        if np.count_nonzero(good) >= 3 and span >= 0.7:
            eligible[i, good] = row[good]
            centered[i] = derivative[i] - np.average(derivative[i, good], weights=row[good])
            spans.append(span)
    support = float(eligible.sum())
    rms = math.sqrt(float(np.sum(eligible * centered ** 2)) / support) if support else 0.0
    total = math.sqrt(float(np.sum(eligible * derivative ** 2)) / support) if support else 0.0
    fraction = rms / max(total, 1e-12)
    identified = rms >= 1.0 and fraction >= 0.05
    confidence = min(1, support / 3) * rms ** 2 / (1 + rms ** 2) * min(1, fraction / 0.1) if identified else 0.0
    return eligible, {"identified": identified, "observations": int(np.count_nonzero(eligible)),
                      "time_span_seconds": max(spans, default=0.0),
                      "centered_derivative_rms_db_per_log_scale": rms,
                      "centered_derivative_fraction": fraction,
                      "conditioning_ratio": total / rms if rms > 1e-12 else None,
                      "confidence": confidence}


def fit_probe_note(triples, reference_offset_db, previous_gain, previous_decay, probe_decay):
    """Fit (velocity, model, probe, reference) Measurements using a local response.

    Loss estimates are deliberately unused: stable beating/rises can carry an
    informative derivative. Floors may constrain gain, never identify decay.
    """
    velocities = [v for v, _, _, _ in triples]
    for _, model, probe, reference in triples:
        validate_measurement_windows((model, probe, reference))
    if len({len(model.times) for _, model, _, _ in triples}) != 1:
        raise ValueError("probe measurements have inconsistent time/window shapes")
    velocity_design = velocity_weights(velocities)
    times = np.array([m.times for _, m, _, _ in triples])
    design = np.repeat(velocity_design, times.shape[1], axis=0)
    gain, decay = np.array(previous_gain, dtype=float), np.array(previous_decay, dtype=float)
    reports = []
    for p in range(PARTIALS):
        derivative, target, tonal_weight = (np.zeros_like(times) for _ in range(3))
        response_weight = np.zeros_like(times)
        log_probe = math.log(probe_decay[p] / previous_decay[p])
        for i, (_, model, probe, reference) in enumerate(triples):
            valid = (model.power[:, p] > 0) & (probe.power[:, p] > 0)
            valid &= (model.confidence[:, p] > 0.2) & (probe.confidence[:, p] > 0.2)
            if abs(log_probe) >= 1e-4:
                derivative[i, valid] = 10 * (np.log10(probe.power[valid, p])
                                            - np.log10(model.power[valid, p])) / log_probe
                response_weight[i, valid] = model.confidence[valid, p] * probe.confidence[valid, p]
            tonal = (reference.confidence[:, p] > 0.2) & (model.confidence[:, p] > 0)
            tonal &= (model.power[:, p] > 0) & (reference.power[:, p] > 0)
            target[i, tonal] = 10 * (np.log10(reference.power[tonal, p])
                                    - np.log10(model.power[tonal, p])) + reference_offset_db
            tonal_weight[i, tonal] = model.confidence[tonal, p] * reference.confidence[tonal, p]
        triple_weight = response_weight * np.array([r.confidence[:, p] for _, _, _, r in triples])
        triple_weight *= tonal_weight > 0
        eligible, condition = response_condition(derivative, triple_weight, times)
        joint = np.zeros(4)
        if condition["identified"]:
            joint_design = np.column_stack((design, derivative.ravel()))
            support = design.T @ eligible.ravel()
            prior = np.r_[0.02 + 0.5 * np.maximum(0, 1 - support), 1.0]
            joint, robust_weight = response_lstsq(joint_design, target.ravel(), eligible.ravel(), prior)
            _, condition = response_condition(derivative, robust_weight.reshape(times.shape), times)
        requested_log = float(joint[3]) if condition["identified"] else 0.0
        bounded_log = float(np.clip(requested_log, *np.log(PROBE_DECAY_STEP)))
        decay[p] = np.clip(previous_decay[p] * math.exp(bounded_log), *DECAY_LIMITS)
        applied_log = math.log(decay[p] / previous_decay[p])
        compensation = derivative * applied_log
        # With a decay update, an unmeasured response cannot be compensated.
        gain_weight = tonal_weight.copy()
        if abs(applied_log) > 1e-12:
            gain_weight *= response_weight > 0
        corrected = target - compensation
        cut_values, cut_weights = [], []
        for i, (_, model, _, reference) in enumerate(triples):
            values, weights = [], []
            for ti in (0, 1):
                if abs(applied_log) > 1e-12 and response_weight[i, ti] <= 0:
                    continue
                upper = reference.power_upper_bound[ti, p]
                if (reference.confidence[ti, p] <= 0.2 and model.confidence[ti, p] > 0.25
                        and reference.location_confidence[p] >= 0.25 and upper > 0 and model.power[ti, p] > 0):
                    ceiling = 10 * math.log10(upper / model.power[ti, p]) + reference_offset_db
                    if ceiling < 0 and ceiling - compensation[i, ti] < 0:
                        values.append(ceiling - compensation[i, ti])
                        weights.append(0.25 * model.confidence[ti, p] * reference.location_confidence[p])
            cut_values.append(max(values) if values else 0.0)
            cut_weights.append(float(np.mean(weights)) if weights else 0.0)
        cut_values, cut_weights = np.array(cut_values), np.array(cut_weights)
        gain_design = np.vstack((design, velocity_design))
        gain_target = np.r_[corrected.ravel(), cut_values]
        weights = np.r_[gain_weight.ravel(), cut_weights]
        support = gain_design.T @ weights
        velocity_fit, _ = response_lstsq(gain_design, gain_target, weights,
                                         0.02 + 0.5 * np.maximum(0, 1 - support))
        by_velocity = gain_weight.mean(axis=1)
        boost_guard = np.clip([by_velocity[velocities.index(v)] / 0.25 if v in velocities else 0 for v in KNOTS], 0, 1)
        raw_update = np.where(velocity_fit > 0, velocity_fit * boost_guard, velocity_fit)
        cut_knots = velocity_design.T @ cut_weights > 0
        raw_update[cut_knots] = np.minimum(raw_update[cut_knots], 0)
        gain[:, p] = np.clip(previous_gain[:, p] + np.clip(raw_update, -GAIN_STEP, GAIN_STEP), *GAIN_LIMITS)
        applied_gain = gain[:, p] - previous_gain[:, p]
        predicted = (design @ applied_gain).reshape(times.shape) + compensation
        residual = target - predicted
        usable = gain_weight > 0
        rms = lambda values: math.sqrt(float(np.sum(gain_weight * values ** 2) / gain_weight.sum())) if np.any(usable) else None
        reports.append({
            "partial": p + 1, "velocities": velocities,
            "gain_confidence_by_velocity": by_velocity.tolist(),
            "residual_db_by_velocity": [float(np.average(row, weights=w)) if np.any(w) else 0.0
                                        for row, w in zip(corrected, gain_weight)],
            "cut_only_upper_db_by_velocity": cut_values.tolist(),
            "cut_only_confidence_by_velocity": cut_weights.tolist(),
            "cut_only_residual_db_before": np.maximum(0, -cut_values).tolist(),
            "cut_only_residual_db_after": np.where(cut_weights > 0,
                np.maximum(0, velocity_design @ applied_gain - cut_values), 0).tolist(),
            "knot_cut_only": cut_knots.tolist(), "knot_boost_evidence_factors": boost_guard.tolist(),
            "velocity_fit_update_db": velocity_fit.tolist(),
            "unclamped_gain_update_db": raw_update.tolist(), "applied_gain_update_db": applied_gain.tolist(),
            "gain_clamped": bool(np.any(abs(raw_update - applied_gain) > 1e-8)),
            "decay_confidence": condition["confidence"],
            "requested_log_decay_update": requested_log, "applied_log_decay_update": applied_log,
            "unclamped_decay_ratio": math.exp(requested_log) if requested_log <= math.log(sys.float_info.max) else None,
            "applied_decay_ratio": float(decay[p] / previous_decay[p]),
            "decay_clamped": abs(requested_log - applied_log) > 1e-8,
            "probe_log_scale_delta": log_probe, "joint_gain_update_db": joint[:3].tolist(),
            "response_condition": condition,
            "derivative_observations": int(np.count_nonzero(response_weight)),
            "derivative_db_per_log_scale": [[float(d) if w > 0 else None for d, w in zip(row, weights)]
                                            for row, weights in zip(derivative, response_weight)],
            "derivative_confidence": response_weight.tolist(),
            "decay_observation_confidence": triple_weight.tolist(),
            "predicted_residual_db": [[float(r) if ok else None for r, ok in zip(row, good)]
                                      for row, good in zip(residual, usable)],
            "predicted_residual_rms_db_before": rms(target), "predicted_residual_rms_db_after": rms(residual),
        })
    return gain, decay, reports


@dataclass(frozen=True)
class Region:
    key: int
    lovel: int
    hivel: int
    sample: Path
    layer: int


def sfz_key(value):
    if re.fullmatch(r"-?\d+", value):
        return int(value)
    match = re.fullmatch(r"([A-Ga-g])([#b]?)(-?\d+)", value)
    if not match:
        raise ValueError(f"invalid SFZ key {value!r}")
    letter, accidental, octave = match.groups()
    return (int(octave) + 1) * 12 + {"C": 0, "D": 2, "E": 4, "F": 5, "G": 7, "A": 9, "B": 11}[letter.upper()] + {"": 0, "#": 1, "b": -1}[accidental]


def parse_sfz(corpus):
    corpus = Path(corpus).resolve()
    sfz = corpus / "SalamanderGrandPianoV3.sfz"
    source = sfz.read_text(encoding="utf-8-sig")
    source = re.sub(r"/\*.*?\*/|//[^\n]*", "", source, flags=re.S)
    if re.search(r"^\s*#", source, re.M):
        raise ValueError(f"{sfz}: SFZ preprocessor directives unsupported")
    global_values, group, regions = {}, {}, []
    for match in re.finditer(r"<([^>]+)>([^<]*)", source):
        section, body = match.groups()
        values = {m.group(1): m.group(2).strip().strip('"') for m in re.finditer(r"(\w+)\s*=\s*(.*?)(?=\s+\w+\s*=|$)", body, re.S)}
        if section == "global":
            global_values = values
            group = {}
        elif section == "group":
            group = values
        elif section == "region":
            values = {**global_values, **group, **values}
            name = values.get("sample", "").replace("\\", "/")
            native = re.fullmatch(r"([A-Ga-g][#b]?-?\d+)v(\d+)\.wav", Path(name).name, re.I)
            if not native or values.get("trigger", "attack") not in ("attack", "first"):
                continue
            key = sfz_key(values.get("pitch_keycenter", values.get("key", "60")))
            low = sfz_key(values.get("lokey", values.get("key", "0")))
            high = sfz_key(values.get("hikey", values.get("key", "127")))
            if key != sfz_key(native[1]) or not low <= key <= high:
                raise ValueError(f"{sfz}: non-native keycenter for {name}")
            if float(values.get("tune", 0)) or float(values.get("transpose", 0)):
                raise ValueError(f"{sfz}: retuned region unsupported: {name}")
            sample = (sfz.parent / name).resolve()
            if not sample.is_relative_to(corpus):
                raise ValueError(f"{sfz}: sample escapes corpus: {name}")
            region = Region(key, int(values.get("lovel", 1)), int(values.get("hivel", 127)), sample, int(native[2]))
            if not (0 <= key <= 127 and 1 <= region.lovel <= region.hivel <= 127 and 1 <= region.layer <= 16):
                raise ValueError(f"{sfz}: invalid region: {name}")
            regions.append(region)
    if not regions:
        raise ValueError(f"{sfz}: no native attack regions")
    return sfz, regions


def select_region(regions, key, velocity):
    matches = [r for r in regions if r.key == key and r.lovel <= velocity <= r.hivel]
    if len(matches) != 1:
        raise ValueError(f"native key {key}, velocity {velocity}: expected one SFZ region, found {len(matches)}")
    return matches[0]


def sha256(path):
    digest = hashlib.sha256()
    with Path(path).open("rb") as source:
        for block in iter(lambda: source.read(1024 * 1024), b""):
            digest.update(block)
    return digest.hexdigest()


def read_previous(path):
    return read_calibration(Path(path).read_text(), str(path))


def read_calibration(text, label):
    notes, seen = {}, set()
    with io.StringIO(text, newline="") as source:
        reader = csv.DictReader(source)
        if reader.fieldnames != list(CSV_FIELDS):
            raise ValueError(f"{label}: expected CSV columns {','.join(CSV_FIELDS)}")
        for row in reader:
            if None in row or any(value is None for value in row.values()):
                raise ValueError(f"{label}: expected six CSV fields per row")
            key, partial = int(row["key"]), int(row["partial"])
            values = np.array([float(row[field]) for field in CSV_FIELDS[2:]])
            if not (0 <= key <= 127 and 1 <= partial <= PARTIALS) or (key, partial) in seen:
                raise ValueError(f"{label}: invalid/duplicate key or partial")
            if not np.isfinite(values).all() or np.any(values[:3] < GAIN_LIMITS[0]) or np.any(values[:3] > GAIN_LIMITS[1]) or not DECAY_LIMITS[0] <= values[3] <= DECAY_LIMITS[1]:
                raise ValueError(f"{label}: non-finite/out-of-range calibration")
            seen.add((key, partial))
            gains, decay = notes.setdefault(key, (np.zeros((3, PARTIALS)), np.ones(PARTIALS)))
            gains[:, partial - 1], decay[partial - 1] = values[:3], values[3]
    if not notes or any(sum(k == key for k, _ in seen) != PARTIALS for key in notes):
        raise ValueError(f"{label}: each key must contain exactly {PARTIALS} partials (legacy 64-row tables are unsupported)")
    return notes


def validate_probe(renders, baseline, probe_dir, probe, previous, factor):
    """Require matching renderer settings and complete embedded calibration tables."""
    if not math.isfinite(factor) or factor <= 0 or factor == 1:
        raise ValueError("--probe-decay-factor must be finite, positive and != 1")
    settings, tables = [], []
    required = {"schema", "mode", "rate_hz", "seconds", "notes", "velocities", "dry", "effects",
                "note_on_sample", "note_off", "block_frames", "design_defaults_plus_overrides"}
    for label, root, info, multiplier in (("baseline", renders, baseline, 1.0),
                                          ("probe", probe_dir, probe, factor)):
        if not isinstance(info, dict) or not required <= info.keys():
            raise ValueError(f"{label}: incomplete renderer settings for decay probe")
        json.dumps(info, allow_nan=False)
        if (info["schema"] != 1 or info["mode"] != "calibration" or info["rate_hz"] != 48000
                or not isinstance(info["seconds"], (int, float)) or info["seconds"] < 4
                or type(info["dry"]) is not bool or not isinstance(info["effects"], str)
                or not isinstance(info["design_defaults_plus_overrides"], dict)
                or info["note_on_sample"] != 0 or info["note_off"] is not None):
            raise ValueError(f"{label}: expected held 48 kHz calibration renders with explicit effects")
        calibration = info.get("calibration")
        if not isinstance(calibration, dict) or not isinstance(calibration.get("csv"), str):
            raise ValueError(f"{label}: render.json must embed calibration.csv text")
        table = read_calibration(calibration["csv"], label)
        if table.keys() != previous.keys():
            raise ValueError(f"{label}: calibration keys differ from --previous")
        for key, (gain, decay) in table.items():
            expected_gain, previous_decay = previous[key]
            # Six-place CSV rounding plus f32 parsing/multiplication precision.
            expected_decay = np.clip(previous_decay * multiplier, *DECAY_LIMITS)
            if (not np.allclose(gain, expected_gain, rtol=2e-7, atol=5.1e-7)
                    or not np.allclose(decay, expected_decay, rtol=2e-7, atol=5.1e-7)):
                raise ValueError(f"{label}: calibration does not match expected gains/decay for key {key}")
        comparable = {k: v for k, v in info.items() if k != "calibration"}
        # The source path and CSV necessarily differ; any other voicing metadata must match.
        comparable["calibration"] = {k: v for k, v in calibration.items() if k not in ("path", "csv")}
        for field in ("notes", "velocities"):
            values = info[field]
            if (not isinstance(values, list) or not values
                    or any(type(v) is not int or not 1 <= v <= 127 for v in values)
                    or len(set(values)) != len(values)):
                raise ValueError(f"{label}: invalid manifest {field}")
            comparable[field] = sorted(values)
        expected_files = {f"note_{key:03}_vel_{velocity:03}.wav"
                          for key in info["notes"] for velocity in info["velocities"]}
        actual_files = {p.name for p in root.glob("note_*_vel_*.wav") if p.is_file()}
        if actual_files != expected_files:
            raise ValueError(f"{label}: WAV notes/velocities do not match render.json")
        settings.append(comparable)
        tables.append(table)
    if settings[0] != settings[1]:
        raise ValueError("baseline/probe renderer settings differ (notes, velocities, voicing or effects)")
    return tables[1]


def write_outputs(out, notes, metadata, summary):
    out = Path(out)
    for key, (gain, decay) in notes.items():
        if not 0 <= key <= 127 or gain.shape != (3, PARTIALS) or decay.shape != (PARTIALS,):
            raise ValueError("invalid calibration shape")
        if not np.isfinite(gain).all() or not np.isfinite(decay).all() or np.any(gain < GAIN_LIMITS[0]) or np.any(gain > GAIN_LIMITS[1]) or np.any(decay < DECAY_LIMITS[0]) or np.any(decay > DECAY_LIMITS[1]):
            raise ValueError("invalid calibration values")
    out.mkdir(parents=True, exist_ok=True)
    with tempfile.TemporaryDirectory(prefix=".fit-", dir=out) as staging:
        staging = Path(staging)
        with (staging / OUTPUTS[0]).open("w", newline="") as destination:
            writer = csv.writer(destination, lineterminator="\n")
            writer.writerow(CSV_FIELDS)
            for key, (gain, decay) in sorted(notes.items()):
                for p in range(PARTIALS):
                    writer.writerow([key, p + 1, *(f"{v:.6f}" for v in gain[:, p]), f"{decay[p]:.6f}"])
        rust = ["// Generated empirical modal voicing; closed-loop validation required.",
                "// Reference: Salamander Grand Piano V3, Alexander Holm, CC BY 3.0.",
                "// Source SHA256 manifest and settings: metadata.json.",
                "use super::CalibrationNote;", "", "pub const DEFAULT_CALIBRATION: &[CalibrationNote] = &["]
        for key, (gain, decay) in sorted(notes.items()):
            rust.extend(["    CalibrationNote {", f"        key: {key},", "        gain_db: ["])
            rust.extend("            [" + ", ".join(f"{v:.6f}" for v in row) + "]," for row in gain)
            rust.extend(["        ],", "        decay_scale: [" + ", ".join(f"{v:.6f}" for v in decay) + "],", "    },"])
        rust.append("];\n")
        (staging / OUTPUTS[1]).write_text("\n".join(rust))
        metadata = {**metadata, "generated_sha256": {name: sha256(staging / name) for name in OUTPUTS[:2]}}
        for name, value in ((OUTPUTS[2], metadata), (OUTPUTS[3], summary)):
            (staging / name).write_text(json.dumps(value, indent=2, sort_keys=True, allow_nan=False) + "\n")
        for name in OUTPUTS:
            os.replace(staging / name, out / name)


def midi_list(value):
    try:
        values = sorted(set(int(v) for v in value.split(",")))
    except ValueError as error:
        raise argparse.ArgumentTypeError("use comma-separated MIDI integers") from error
    if not values or min(values) < 1 or max(values) > 127:
        raise argparse.ArgumentTypeError("MIDI values must be in 1..127")
    return values


def run(args):
    renders, corpus, out = (Path(p).resolve() for p in (args.renders, args.corpus, args.out))
    probe_dir = Path(args.decay_probe).resolve() if args.decay_probe else None
    if (probe_dir is None) != (args.probe_decay_factor is None) or (probe_dir and not args.previous):
        raise ValueError("--decay-probe and --probe-decay-factor require each other and --previous")
    if out.is_relative_to(renders) or out.is_relative_to(corpus) or (probe_dir and out.is_relative_to(probe_dir)):
        raise ValueError("--out must be outside renders, probe and corpus to preserve sources")
    previous_path = Path(args.previous).resolve() if args.previous else None
    if previous_path and previous_path in [(out / name).resolve() for name in OUTPUTS]:
        raise ValueError("--out must not overwrite --previous")
    render_json = renders / "render.json"
    render_info = json.loads(render_json.read_text())
    if not isinstance(render_info, dict):
        raise ValueError("render.json must contain a JSON object")
    sfz, regions = parse_sfz(corpus)
    native_keys = sorted({region.key for region in regions})
    keys = args.notes or native_keys
    if not args.notes and len(native_keys) != 30:
        raise ValueError(f"expected 30 native keys, found {len(native_keys)}; use --notes for a pilot")
    if any(key not in native_keys for key in keys):
        raise ValueError("--notes contains a non-native SFZ keycenter")
    expected = sorted(set(KNOTS) | set(args.velocities))
    inputs = {}
    for key in keys:
        velocities = set(expected)
        for path in renders.glob(f"note_{key:03}_vel_*.wav"):
            match = re.fullmatch(r"note_\d{3}_vel_(\d{3})\.wav", path.name)
            if not match or not 1 <= int(match[1]) <= 127:
                raise ValueError(f"invalid render filename: {path.name}")
            velocities.add(int(match[1]))
        for velocity in sorted(velocities):
            region = select_region(regions, key, velocity)
            inputs[key, velocity] = (renders / f"note_{key:03}_vel_{velocity:03}.wav", region)
    previous = read_previous(previous_path) if previous_path else {}
    if previous and any(key not in previous for key in keys):
        raise ValueError("--previous must contain every fitted key")
    probe_info, probe_table = None, None
    if probe_dir:
        probe_info = json.loads((probe_dir / "render.json").read_text())
        probe_table = validate_probe(renders, render_info, probe_dir, probe_info, previous, args.probe_decay_factor)
    anchor = {"method": "fixed", "offset_db": args.reference_offset_db}
    if previous_path and (args.reference_offset_db is None or (probe_dir and (previous_path.parent / "metadata.json").is_file())):
        previous_metadata = previous_path.parent / "metadata.json"
        prior = json.loads(previous_metadata.read_text())
        if prior.get("generated_sha256", {}).get("calibration.csv") != sha256(previous_path):
            if probe_dir:
                raise ValueError("previous metadata does not identify this CSV; decay probe requires the original anchor provenance")
            raise ValueError("previous metadata does not identify this CSV; supply --reference-offset-db explicitly")
        if probe_dir and args.reference_offset_db is not None and not math.isclose(args.reference_offset_db, prior["anchor"]["offset_db"], abs_tol=1e-9, rel_tol=0):
            raise ValueError("decay probe must reuse the previous global reference offset")
        anchor = {**prior["anchor"], "reused_from_previous": True}
    elif args.reference_offset_db is None:
        region = select_region(regions, 60, 68)
        inputs.setdefault((60, 68), (renders / "note_060_vel_068.wav", region))
        anchor = {"method": "C4_68_stereo_rms_0.05_to_0.45s", "key": 60, "velocity": 68, "offset_db": None}
    if (probe_dir and anchor["offset_db"] is None) or (anchor["offset_db"] is not None and not math.isfinite(anchor["offset_db"])):
        raise ValueError("reference offset must be finite")
    missing = [str(path) for model, region in inputs.values() for path in (model, region.sample) if not path.is_file()]
    if missing:
        raise ValueError("missing required inputs (no pitch/layer substitution):\n" + "\n".join(sorted(set(missing))))
    source_files = {sfz, *(region.sample for _, region in inputs.values())}
    if (corpus / "README").is_file():
        source_files.add(corpus / "README")
    manifest = lambda paths, root: [{"path": str(path.relative_to(root)), "sha256": sha256(path), "bytes": path.stat().st_size} for path in sorted(paths)]
    source_manifest = manifest(source_files, corpus)
    render_manifest = manifest({render_json, *(model for model, _ in inputs.values())}, renders)
    probe_manifest = manifest({probe_dir / "render.json", *(probe_dir / model.name for model, _ in inputs.values())}, probe_dir) if probe_dir else None

    def measure(path, key, model):
        samples, rate, fmt = read_wav(path)
        if rate != 48000 or fmt["channels"] != 2 or fmt["encoding"] != (3 if model else 1) or fmt["bits"] != (32 if model else 24):
            raise ValueError(f"{path}: expected stereo 48 kHz {'float32 render' if model else 'PCM24 reference'}")
        if model and len(samples) < rate * 4 - 1:
            raise ValueError(f"{path}: expected held render of at least 4 seconds")
        return analyze(samples, rate, key)

    cached_anchor = None
    if anchor["offset_db"] is None:
        model_path, region = inputs[60, 68]
        model, reference = measure(model_path, 60, True), measure(region.sample, 60, False)
        anchor["offset_db"] = 20 * math.log10(model.report["early_stereo_rms"] / reference.report["early_stereo_rms"])
        anchor["model_rms"], anchor["reference_rms"] = model.report["early_stereo_rms"], reference.report["early_stereo_rms"]
        anchor["original_model_sha256"] = sha256(model_path)
        anchor["original_reference_sha256"] = sha256(region.sample)
        cached_anchor = (model, reference)
    notes, reports = {}, []
    for key in keys:
        pairs, triples, observations, reference_cache = [], [], [], {}
        for (input_key, velocity), (model_path, region) in sorted(inputs.items()):
            if input_key != key:
                continue
            if (key, velocity) == (60, 68) and cached_anchor is not None:
                model, reference = cached_anchor
            else:
                model = measure(model_path, key, True)
                if region.sample not in reference_cache:
                    reference_cache[region.sample] = measure(region.sample, key, False)
                reference = reference_cache[region.sample]
            pairs.append((velocity, model, reference))
            observations.append({"velocity": velocity, "layer": region.layer,
                                 "reference": str(region.sample.relative_to(corpus)),
                                 "model": model.report, "reference_measurement": reference.report})
            if probe_dir:
                probe = measure(probe_dir / model_path.name, key, True)
                triples.append((velocity, model, probe, reference))
                observations[-1]["probe_measurement"] = probe.report
        if probe_dir:
            gain, decay, partial_report = fit_probe_note(triples, anchor["offset_db"], *previous[key], probe_table[key][1])
        else:
            gain, decay, partial_report = fit_note(pairs, anchor["offset_db"], *previous.get(key, (None, None)))
        notes[key] = gain, decay
        reports.append({"key": key, "inputs": observations, "partials": partial_report})
        print(f"fitted native key {key}: {len(pairs)} velocities", file=sys.stderr)
    metadata = {
        "schema_version": 1, "algorithm": "empirical-stereo-modal-voicing-v1",
        "script_sha256": sha256(__file__), "numpy_version": np.__version__, "python_version": sys.version,
        "reference_attribution": {"title": "Salamander Grand Piano V3", "author": "Alexander Holm", "license": "CC BY 3.0", "license_url": "https://creativecommons.org/licenses/by/3.0/"},
        "source_manifest": source_manifest, "render_manifest": render_manifest,
        "render_metadata": render_info, "anchor": anchor, "keys": keys,
        "expected_velocities": expected, "velocity_knots": KNOTS, "partials": PARTIALS,
        "windows_seconds": TIMES.tolist(), "window_durations_seconds": WINDOW_SECONDS.tolist(),
        "gain_limits_db": GAIN_LIMITS, "gain_step_limit_db": GAIN_STEP,
        "decay_limits": DECAY_LIMITS, "decay_step_limits": DECAY_STEP,
        "measurement_settings": {
            "onset_block_seconds": 0.001, "onset_peak_search_seconds": 0.5,
            "onset_relative_threshold_db": -40, "identification_seconds": [0.1, 1.2],
            "early_identification_seconds": [0, EARLY_SECONDS],
            "early_match_radius": "min(0.23 * harmonic gap, 2.5 * max(2 / early duration, 0.003 * predicted Hz))",
            "f0_search_cents": [-45, 45], "B_search": [0, 0.01],
            "line_prominence_min_db": 15, "band_snr_fade_db": [10, 25],
            "relative_band_power_fade_db": [-60, -40],
            "band_half_width": "min(0.23 * harmonic gap, max(2.5 / window duration, 0.008 * center Hz))",
            "gain_velocity_zero_prior": "0.02 + 0.5 * max(0, 1 - knot support)",
            "positive_gain_knot_evidence_full_at": 0.25,
            "power_upper_bound": "predicted harmonic band total power plus one local-floor allowance",
            "cut_only_weight_scale": 0.25,
            "cut_only_location_confidence_min": 0.25,
            "cut_only_model_confidence_min": 0.25,
            "cut_only_reference_confidence_max": 0.2,
            "decay_neighborhood_strength_max": 0.15,
        },
        "previous_sha256": sha256(previous_path) if previous_path else None,
        "reference_target": "native recorded PCM; SFZ used for key/layer selection, no SFZ gain or retuning",
        "acceptance": "unvalidated; requires held-out and closed-loop rerender measurement",
    }
    if probe_dir:
        metadata.update({"algorithm": "empirical-stereo-modal-voicing-probe-v1",
                         "probe_manifest": probe_manifest, "probe_metadata": probe_info,
                         "probe_decay_factor": args.probe_decay_factor, "decay_step_limits": PROBE_DECAY_STEP})
        metadata["measurement_settings"].update({
            "response": "(probe_power_db - model_power_db) / log(actual_probe_scale / previous_scale)",
            "probe_log_scale_delta_min": 1e-4, "response_tonal_confidence_min": 0.2,
            "response_time_support_weight_min": 0.05,
            "response_min_windows_per_velocity": 3, "response_min_time_span_seconds": 0.7,
            "response_centered_rms_min_db_per_log_scale": 1.0, "response_centered_fraction_min": 0.05,
            "response_huber_db": 3.0, "response_irls_iterations": 12, "response_log_decay_zero_prior": 1.0,
            "decay_neighborhood_strength_max": 0.0,
        })
    all_partials = [p for report in reports for p in report["partials"]]
    summary = {"acceptance": metadata["acceptance"], "anchor": anchor,
               "fitted_keys": len(notes), "fitted_inputs": sum(len(r["inputs"]) for r in reports),
               "gain_clamped_partials": sum(p["gain_clamped"] for p in all_partials),
               "decay_clamped_partials": sum(p["decay_clamped"] for p in all_partials),
               "unsupported_partials": sum(not any(p["gain_confidence_by_velocity"] + p["cut_only_confidence_by_velocity"]) for p in all_partials),
               "notes": reports}
    write_outputs(out, notes, metadata, summary)
    print(f"Wrote {len(notes)} keys / {len(notes) * PARTIALS} rows to {out}; closed-loop validation required.")


def main(argv=None):
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--renders", required=True, type=Path)
    parser.add_argument("--corpus", required=True, type=Path)
    parser.add_argument("--out", required=True, type=Path, help="explicit output directory; only four named outputs may be replaced")
    parser.add_argument("--previous", type=Path)
    parser.add_argument("--decay-probe", type=Path, help="matched renders with previous decay scales perturbed; requires --previous and --probe-decay-factor")
    parser.add_argument("--probe-decay-factor", type=float, help="positive non-unit multiplier used for every probe decay scale before clamping (typically 0.7)")
    parser.add_argument("--velocities", type=midi_list, default=list(KNOTS), help="comma-separated required velocities; knots are always required; available extra layers also contribute")
    parser.add_argument("--notes", type=midi_list, help="comma-separated native MIDI keys for a pilot; default all 30")
    parser.add_argument("--reference-offset-db", type=float, help="fixed shared dB offset added to reference; otherwise anchor C4/68, or reuse previous metadata")
    args = parser.parse_args(argv)
    try:
        run(args)
    except (OSError, ValueError, KeyError, struct.error) as error:
        parser.exit(2, f"fit_voicing: {error}\n")


if __name__ == "__main__":
    main()