secs-engine 0.10.0

//! Per-channel room correction, in upstream arithmetic and decision order.
//! Upstream: `_correction`, its masks, and `compute`'s per-main loop.

use crate::{
    C, EngineConfig, EngineError, Tuning, decay,
    fit::{self, FitRequest, Fitter},
    grid::{BINS, N_FFT, SR_F, design_freqs},
    input::{Peq, PhasePositions},
    levels, measurement,
    notes::Note,
    setup::Setup,
};
use secs_dsp::{
    complex_div as divide, fft, interp, lowrate,
    minphase::minimum_phase,
    pairwise_sum,
    response::{self, magnitude},
    smooth::{self, Fraction},
};
use std::{f64::consts::PI, sync::Arc};

pub struct RcInput<'a> {
    pub role: &'a str,
    pub comb: &'a [(u32, Arc<[C]>)],
    pub main_position: u32,
    pub level_db: f64,
    pub setup: &'a Setup,
    pub config: &'a EngineConfig,
    pub tuning: &'a Tuning,
    pub positions: &'a [u32],
}
#[derive(Debug, Clone)]
pub struct RcOutput {
    pub fir: Vec<f64>,
    pub response: Vec<C>,
    pub peq: Vec<Peq>,
    pub phase_strength: f64,
    pub phase_band: Option<f64>,
    pub design_db: Vec<f64>,
    pub notes: Vec<Note>,
}
pub struct ChannelRc {
    pub out: RcOutput,
    pub phase_frac_used: Option<f64>,
    pub gd_cuts: Vec<Peq>,
    pub response: Vec<C>,
}

// The trace observes the same implementation; it never supplies decisions.
#[derive(Default, Debug)]
pub(crate) struct Trace {
    pub no_boost: Vec<bool>,
    pub long_mode: Vec<bool>,
    pub peaks: Vec<(u32, usize)>,
    pub ladder: Vec<(f64, f64, f64)>,
}
fn main_response<'a>(input: &'a RcInput<'_>) -> &'a [C] {
    &input
        .comb
        .iter()
        .find(|(p, _)| *p == input.main_position)
        .expect("measured main position")
        .1
}
fn smooth(x: &[f64], fraction: f64) -> Vec<f64> {
    smooth::smooth_power(design_freqs(), x, Fraction::Scalar(fraction)).expect("design grid")
}
/// The magnitude estimate used both by RC fitting and the agreement probe.
fn smoothed_db(power: &[f64]) -> Vec<f64> {
    let f = design_freqs();
    let frac = secs_dsp::grid::frac_for(f);
    let sm = smooth::smooth_power(f, power, Fraction::PerBin(&frac)).expect("design grid");
    let sm3 = smooth(power, 3.0);
    sm.iter()
        .zip(sm3)
        .map(|(sm, sm3)| (10.0 * (sm + 1e-30).log10()).max(10.0 * (sm3 + 1e-30).log10() - 3.0))
        .collect()
}

/// Evaluate the actual pass-1 RC response against each position's target error.
/// Position weights apply to gains, not to spectra averaged before scoring.
pub(crate) fn agreement_limit(input: &RcInput<'_>, correction: &[C]) -> f64 {
    let f = design_freqs();
    let width = 2.0_f64.powf(1.0 / 6.0);
    let bands: Vec<_> = levels::RC_BAND_CENTRES
        .iter()
        .map(|fc| f.partition_point(|f| *f < fc / width)..f.partition_point(|f| *f <= fc * width))
        .collect();
    let target = levels::target_db(
        &input.config.target,
        &f.iter().map(|v| v.max(1.0)).collect::<Vec<_>>(),
    );
    // Use the same power-to-dB floor as the measured curve. The fitting
    // helper `db` adds 1e-12 to magnitude, making unity look like a tiny
    // improvement for a negative error instead of the required zero gain.
    let correction: Vec<_> = correction
        .iter()
        .map(|&h| {
            let m = magnitude(h);
            10.0 * (m * m + 1e-30).log10()
        })
        .collect();
    let weights = levels::weights(input.positions, input.main_position);
    let mut gains = [0.0; 14];
    let mut total = 0.0;
    for (position, h) in input.comb {
        let weight = weights[input
            .positions
            .iter()
            .position(|p| p == position)
            .expect("engine position")];
        let power: Vec<_> = h
            .iter()
            .map(|v| {
                let m = magnitude(*v);
                m * m
            })
            .collect();
        let measured = smoothed_db(&power);
        for (gain, band) in gains.iter_mut().zip(&bands) {
            let mut before = 0.0;
            let mut after = 0.0;
            for i in band.clone() {
                let error = measured[i] - target[i] - input.level_db;
                before += error;
                after += error + correction[i];
            }
            *gain += weight * (before.abs() - after.abs()) / band.len() as f64;
        }
        total += weight;
    }
    let mut limit = 100.0;
    for (fc, gain) in levels::RC_BAND_CENTRES.iter().zip(gains) {
        if (gain / total).partial_cmp(&0.0) != Some(std::cmp::Ordering::Greater) {
            break;
        }
        limit = fc * width;
    }
    // Clamp only after taking the median over mains. Clamping each main here
    // would lower the median of an even-sized set near the 2 kHz boundary.
    limit
}

/// Sorted union of the resolved primary band and the manually set extra bands.
/// Keep legacy primary bounds untouched when there are no extras.
pub(crate) fn band_set(lo: f64, hi: f64, extra: &[(f64, f64)]) -> Vec<(f64, f64)> {
    let mut bands = vec![(lo, hi)];
    bands.extend_from_slice(extra);
    merge(bands)
}

/// Sort by lower edge and merge intervals that overlap or touch.
fn merge(mut bands: Vec<(f64, f64)>) -> Vec<(f64, f64)> {
    bands.sort_by(|a, b| a.0.total_cmp(&b.0));
    let mut merged: Vec<(f64, f64)> = Vec::new();
    for (lo, hi) in bands {
        if let Some(last) = merged.last_mut()
            && lo <= last.1
        {
            last.1 = last.1.max(hi);
        } else {
            merged.push((lo, hi));
        }
    }
    merged
}

pub(crate) fn in_bands(f: f64, bands: &[(f64, f64)]) -> bool {
    bands.iter().any(|&(lo, hi)| f >= lo && f <= hi)
}

/// A geometric grid of `n` points per band, from `max(lo, floor)` to `hi`.
/// A lone band is the legacy grid as it was, falling when an automatic upper
/// limit lies below `rc_lo`. With extra bands an interval that holds no
/// frequency (`max(lo, floor) >= hi`) gets no points, so every point lies in
/// a band.
pub(crate) fn band_grid(bands: &[(f64, f64)], floor: f64, n: usize) -> Vec<f64> {
    if let [(lo, hi)] = bands {
        return geomspace(lo.max(floor), *hi, n);
    }
    bands
        .iter()
        .filter(|&&(lo, hi)| lo.max(floor) < hi)
        .flat_map(|&(lo, hi)| geomspace(lo.max(floor), hi, n))
        .collect()
}

/// What a Zero-mode PEQ fit over `lo..hi` weighs (ADR 0024). Without an extra
/// band inside `lo..hi` the fit keeps its legacy weights (`None`), so an extra
/// band elsewhere changes nothing in it. With one, the fit weighs where the
/// magnitude ramps reach: the primary band from `lo` to `rc_hi·√2` and each
/// extra band from `lo/√2` to `hi·√2`; outside their union it seeds nothing
/// and weighs 0.15.
pub(crate) fn fit_reach(setup: &Setup, lo: f64, hi: f64) -> Option<Vec<(f64, f64)>> {
    if !setup
        .rc_bands
        .iter()
        .any(|&(band_lo, band_hi)| band_lo <= hi && band_hi >= lo)
    {
        return None;
    }
    let sqrt2 = 2_f64.sqrt();
    Some(merge(
        std::iter::once((lo, setup.rc_hi * sqrt2))
            .chain(
                setup
                    .rc_bands
                    .iter()
                    .map(|&(band_lo, band_hi)| (band_lo / sqrt2, band_hi * sqrt2)),
            )
            .filter(|&(lo, hi)| lo < hi)
            .collect(),
    ))
}

pub(crate) fn magnitude_mask(f: f64, bands: &[(f64, f64)]) -> f64 {
    let lf = f.max(1e-3).log2();
    bands
        .iter()
        .map(|&(lo, hi)| {
            let hi_edge = (hi * 2_f64.sqrt()).min(SR_F / 2.0 * 0.99);
            ((lf - (lo / 2_f64.sqrt()).log2()) / 0.5).clamp(0.0, 1.0)
                * ((hi_edge.log2() - lf) / (hi_edge / hi).log2().max(1e-3)).clamp(0.0, 1.0)
        })
        .fold(0.0, f64::max)
}

pub(crate) fn phase_fade(f: f64, setup: &Setup, bands: &[(f64, f64)]) -> f64 {
    let fade_hz = setup.rc_hi * 2_f64.sqrt();
    let primary = if f < 15.0 {
        (f / 15.0).clamp(0.0, 1.0)
    } else {
        ((fade_hz.log2() - f.max(1e-3).log2()) / 0.5).clamp(0.0, 1.0)
    };
    if setup.rc_bands.is_empty() {
        return primary;
    }
    let lf = f.max(1e-3).log2();
    bands
        .iter()
        .map(|&(lo, hi)| {
            ((lf - (lo / 2_f64.sqrt()).log2()) / 0.5).clamp(0.0, 1.0)
                * (((hi * 2_f64.sqrt()).log2() - lf) / 0.5).clamp(0.0, 1.0)
        })
        .fold(primary, f64::max)
}

/// Where the realizer measures the FIR's magnitude error when it picks the
/// phase strength: the primary band from `max(20, lo_fir)` to `rc_hi·√2`, as
/// without extra bands, and each extra band from its lower edge to its upper
/// edge ×√2 below Nyquist (ADR 0024).
pub(crate) fn phase_error_band(f: f64, lo_fir: f64, setup: &Setup) -> bool {
    f >= 20_f64.max(lo_fir)
        && (f <= setup.rc_hi * 2_f64.sqrt()
            || setup
                .rc_bands
                .iter()
                .any(|&(lo, hi)| f >= lo && f <= (hi * 2_f64.sqrt()).min(SR_F / 2.0 * 0.99)))
}

fn geomspace(lo: f64, hi: f64, n: usize) -> Vec<f64> {
    let start = lo.log10();
    let step = (hi.log10() - start) / (n - 1) as f64;
    let mut g: Vec<_> = (0..n)
        .map(|i| 10_f64.powf(i as f64 * step + start))
        .collect();
    g[0] = lo;
    g[n - 1] = hi;
    g
}
fn db(h: C) -> f64 {
    20.0 * (magnitude(h) + 1e-12).log10()
}
fn multiply(a: &[C], b: &[C]) -> Vec<C> {
    a.iter().zip(b).map(|(a, b)| a * b).collect()
}
fn rms(v: &[f64]) -> f64 {
    (pairwise_sum(&v.iter().map(|v| v * v).collect::<Vec<_>>()) / v.len() as f64).sqrt()
}

fn no_boost_mask(input: &RcInput<'_>) -> Vec<bool> {
    let f = design_freqs();
    let bands = band_set(input.setup.rc_lo, input.setup.rc_hi, &input.setup.rc_bands);
    let g = band_grid(&bands, 20.0, 1200);
    let gd =
        response::gd_curve(main_response(input), f, &g, Fraction::Scalar(24.0)).expect("GD grid");
    let mut mask = vec![false; BINS];
    let width = 2_f64.powf(1.0 / 12.0);
    for i in 0..g.len() {
        if gd[i] - decay::median(&gd[i.saturating_sub(60)..(i + 60).min(g.len())])
            > input.tuning.noboost_ms
        {
            for (m, &f) in mask.iter_mut().zip(f) {
                *m |= f >= g[i] / width && f <= g[i] * width;
            }
        }
    }
    mask
}
fn long_mode_mask(input: &RcInput<'_>, notes: &mut Vec<Note>) -> Vec<bool> {
    let f = design_freqs();
    let x = lowrate::to_low(main_response(input), f, 300).expect("low-rate grid");
    let bands = band_set(input.setup.rc_lo, input.setup.rc_hi, &input.setup.rc_bands);
    let modes: Vec<_> = decay::late_peaks(&[&x], 8)
        .into_iter()
        .filter(|&p| p <= 160.0 && in_bands(p, &bands))
        .filter_map(|p| {
            decay::mode_t60(&x, 600.0, p, 1.5)
                .filter(|t| *t != 0.0)
                .map(|t| (p, t))
        })
        .collect();
    let mut mask = vec![false; BINS];
    if modes.len() < 2 {
        return mask;
    }
    let threshold = (1.3 * decay::median(&modes.iter().map(|v| v.1).collect::<Vec<_>>())).max(0.3);
    let mut selected: Vec<_> = modes.into_iter().filter(|v| v.1 >= threshold).collect();
    let width = 2_f64.powf(1.0 / 24.0);
    for &(p, _) in &selected {
        for (m, &f) in mask.iter_mut().zip(f) {
            *m |= f >= p / width && f <= p * width;
        }
    }
    if mask.iter().any(|m| *m) {
        selected.sort_by(|a, b| a.0.total_cmp(&b.0));
        notes.push(Note::LongModeBoost {
            cap_db: input.tuning.long_boost_db,
            modes: selected,
        });
    }
    mask
}
#[allow(clippy::too_many_arguments)]
fn fit_request(
    input: &RcInput<'_>,
    g: &[f64],
    c: &[f64],
    n: usize,
    lo: f64,
    hi: f64,
    f_max: Option<f64>,
    fitter: &mut Fitter<'_>,
) -> Vec<Peq> {
    let reach = (input.config.taps == 0)
        .then(|| fit_reach(input.setup, lo, hi))
        .flatten();
    fitter(&FitRequest {
        g,
        corr_db: c,
        n,
        lo,
        hi,
        boost: input.config.max_boost_db,
        ring_ms: input.tuning.boost_ring_ms,
        bands: reach.as_deref(),
        f_max,
    })
}
fn unwrap(h: &[C]) -> Vec<f64> {
    let angle: Vec<_> = h
        .iter()
        .map(|h| {
            let clean = |x: f64| {
                if x.is_nan() {
                    0.0
                } else {
                    x.clamp(-f64::MAX, f64::MAX)
                }
            };
            clean(h.im).atan2(clean(h.re))
        })
        .collect();
    let mut out = angle.clone();
    let mut correction = 0.0;
    for i in 1..out.len() {
        let d = angle[i] - angle[i - 1];
        let mut dm = (d + PI).rem_euclid(2.0 * PI) - PI;
        if dm == -PI && d > 0.0 {
            dm = PI;
        }
        correction += if d.abs() < PI { 0.0 } else { dm - d };
        out[i] += correction;
    }
    out
}
fn gradient(phase: &[f64]) -> Vec<f64> {
    let f = design_freqs();
    let x: Vec<_> = f.iter().map(|v| 2.0 * PI * v.max(f[1])).collect();
    let dx: Vec<_> = x.windows(2).map(|x| x[1] - x[0]).collect();
    let n = phase.len();
    let mut tau = vec![0.0; n];
    tau[0] = -(phase[1] - phase[0]) / dx[0];
    tau[n - 1] = -(phase[n - 1] - phase[n - 2]) / dx[n - 2];
    for i in 1..n - 1 {
        let (a, b) = (dx[i - 1], dx[i]);
        tau[i] = -((-b / (a * (a + b))) * phase[i - 1]
            + ((b - a) / (a * b)) * phase[i]
            + (a / (b * (a + b))) * phase[i + 1]);
    }
    for v in &mut tau {
        if !v.is_finite() {
            *v = 0.0;
        }
    }
    tau
}
fn truncate(spec: &[C], taps: usize) -> (Vec<f64>, Vec<C>) {
    let mut h = fft::irfft(spec, N_FFT);
    h.truncate(taps);
    let k = (taps as f64 * 0.1) as usize;
    for (i, v) in h[taps - k..].iter_mut().enumerate() {
        *v *= 0.5 * (1.0 + (PI * i as f64 / k as f64).cos());
    }
    let f = measurement::spectrum(&h);
    (h, f)
}
struct Realizer<'a> {
    excess: &'a [f64],
    band: &'a [bool],
    lift: f64,
    taps: usize,
}
impl Realizer<'_> {
    fn realize(&self, frac: f64, m: &[C]) -> (Vec<f64>, Vec<C>, f64, f64) {
        let mut sum = 0.0;
        let spec: Vec<_> = m
            .iter()
            .zip(self.excess)
            .map(|(&m, &e)| {
                sum += self.lift - e.clamp(-0.005, self.lift * frac);
                m * C::from_polar(1.0, -2.0 * PI * sum * design_freqs()[1])
            })
            .collect();
        let (h, f) = truncate(&spec, self.taps);
        let err: Vec<_> = f
            .iter()
            .zip(m)
            .zip(self.band)
            .filter(|(_, b)| **b)
            .map(|((&f, &m), _)| db(f) - db(m))
            .collect();
        let rms = rms(&err);
        let max = err.iter().map(|v| v.abs()).fold(0.0, f64::max);
        (h, f, rms, max)
    }
}
pub fn correction(input: &RcInput<'_>, phase_frac: f64, fitter: &mut Fitter<'_>) -> RcOutput {
    correction_traced(input, phase_frac, fitter, &mut Trace::default())
}

/// `Err(Cancelled)` once `cancelled()` reads true (TD-49). RC polls it at the
/// boundaries of its own stages (masks, each PEQ fit, each phase position,
/// each ladder step and trim pass, the second phase candidate and the GD
/// cuts); the checks change no arithmetic or order.
fn check(cancelled: &dyn Fn() -> bool) -> Result<(), EngineError> {
    if cancelled() {
        Err(EngineError::Cancelled)
    } else {
        Ok(())
    }
}
// The first PEQ fit and position-common excess phase do not depend on the
// candidate phase fraction. Each candidate owns only the values it modifies.
struct FirPreparation {
    peq: Vec<Peq>,
    p: Vec<C>,
    m: Vec<C>,
    u: Vec<C>,
    nb: Vec<bool>,
    g: Vec<f64>,
    lo_fir: f64,
    design_db: Vec<f64>,
    notes: Vec<Note>,
}

enum RcPreparation {
    Complete(RcOutput),
    Fir(FirPreparation),
}

pub(crate) fn correction_traced(
    input: &RcInput<'_>,
    phase_frac: f64,
    fitter: &mut Fitter<'_>,
    trace: &mut Trace,
) -> RcOutput {
    // A probe that never reads true cannot return `Cancelled`.
    let never = || false;
    let prepared = prepare_correction(input, fitter, trace, &never).expect("never cancelled");
    finish_correction(input, phase_frac, &prepared, fitter, trace, &never).expect("never cancelled")
}

fn prepare_correction(
    input: &RcInput<'_>,
    fitter: &mut Fitter<'_>,
    trace: &mut Trace,
    cancelled: &dyn Fn() -> bool,
) -> Result<RcPreparation, EngineError> {
    check(cancelled)?;
    let f = design_freqs();
    let (s, c, t) = (input.setup, input.config, input.tuning);
    let weights = levels::freq_weights(input.positions, input.main_position, t.main_weight_lo);
    let mut power = vec![0.0; BINS];
    let mut total = vec![0.0; BINS];
    for (pos, h) in input.comb {
        let w = &weights[input
            .positions
            .iter()
            .position(|p| p == pos)
            .expect("engine position")];
        for i in 0..BINS {
            let m = magnitude(h[i]);
            power[i] += w[i] * (m * m);
            total[i] += w[i];
        }
    }
    for (p, w) in power.iter_mut().zip(total) {
        *p /= w;
    }
    let sm = smoothed_db(&power);
    let target = c
        .target
        .sample(&f.iter().map(|v| v.max(1.0)).collect::<Vec<_>>());
    let mut corr: Vec<_> = (0..BINS)
        .map(|i| {
            let sm = sm[i];
            (target[i] + input.level_db - sm).clamp(
                -25.0,
                if f[i] > 16000.0 {
                    c.max_boost_db.min(2.0)
                } else {
                    c.max_boost_db
                },
            )
        })
        .collect();
    let mut notes = Vec::new();
    let mut nb = no_boost_mask(input);
    let rolloff = input
        .config
        .low_protect
        .then(|| crate::low_protect::rolloff(input.comb, input.main_position, None))
        .flatten();
    if let Some(hz) = rolloff {
        for (masked, &f) in nb.iter_mut().zip(f) {
            *masked |= f < hz;
        }
    }
    let lm = long_mode_mask(input, &mut notes);
    check(cancelled)?;
    trace.no_boost = nb.clone();
    trace.long_mode = lm.clone();
    let bands = band_set(s.rc_lo, s.rc_hi, &s.rc_bands);
    for i in 0..BINS {
        if nb[i] {
            corr[i] = corr[i].min(0.0);
        }
        if lm[i] {
            corr[i] = corr[i].min(t.long_boost_db);
        }
        corr[i] *= magnitude_mask(f[i], &bands);
    }
    let design_db = corr.clone();
    let mut peq = Vec::new();
    let mut p = vec![C::new(1.0, 0.0); BINS];
    if c.taps == 0 {
        let g = geomspace(15.0, 20000.0, 900);
        let cg: Vec<_> = g.iter().map(|&g| interp(g, f, &corr)).collect();
        let low: Vec<_> = g
            .iter()
            .zip(&cg)
            .map(|(&g, &v)| if g <= 500.0 { v } else { 0.0 })
            .collect();
        peq = fit_request(
            input,
            &g,
            &low,
            s.peq_bands / 2,
            15.0_f64.max(s.rc_lo / 1.4),
            400.0,
            None,
            fitter,
        );
        check(cancelled)?;
        let rest: Vec<_> = g
            .iter()
            .zip(&cg)
            .zip(response::peq_response(&g, &peq))
            .map(|((&g, &v), p)| if g >= 300.0 { v - db(p) } else { 0.0 })
            .collect();
        peq.extend(fit_request(
            input,
            &g,
            &rest,
            s.peq_bands - s.peq_bands / 2,
            400.0,
            bands.last().expect("primary band").1.min(20000.0),
            None,
            fitter,
        ));
        check(cancelled)?;
        if let Some(hz) = rolloff {
            crate::low_protect::clamp_centres(&mut peq, hz);
        }
        return Ok(RcPreparation::Complete(RcOutput {
            fir: vec![1.0],
            response: response::peq_response(f, &peq),
            peq,
            phase_strength: 0.0,
            phase_band: None,
            design_db,
            notes,
        }));
    }
    let g = geomspace(15.0, 20000.0, 700);
    let lo_fir = 1.5 * SR_F / c.taps as f64;
    if s.peq_bands > 0 {
        let cg: Vec<_> = g.iter().map(|&g| interp(g, f, &corr)).collect();
        peq = fit_request(
            input,
            &g,
            &cg,
            s.peq_bands,
            20.0,
            500.0,
            Some(s.rc_hi),
            fitter,
        );
        check(cancelled)?;
        if let Some(hz) = rolloff {
            crate::low_protect::clamp_centres(&mut peq, hz);
        }
        p = response::peq_response(f, &peq);
        for i in 0..BINS {
            corr[i] =
                (corr[i] - db(p[i])) * ((f[i] - lo_fir / 2.0) / (lo_fir / 2.0)).clamp(0.0, 1.0);
        }
    }
    let m = minimum_phase(
        &corr
            .iter()
            .map(|v| 10_f64.powf(v / 20.0))
            .collect::<Vec<_>>(),
    )
    .expect("magnitude");
    let weights = levels::weights(input.positions, input.main_position);
    let frac = secs_dsp::grid::frac_for(f);
    let mut u = vec![C::new(0.0, 0.0); BINS];
    let mut total = 0.0;
    for (pos, h) in input.comb {
        if t.phase_positions == PhasePositions::Main && *pos != input.main_position {
            continue;
        }
        check(cancelled)?;
        let w = weights[input.positions.iter().position(|p| p == pos).unwrap()];
        total += w;
        let mut reference = multiply(&multiply(h, &p), &m);
        let ir = fft::irfft(&reference, N_FFT);
        let peak = decay::first_max(&ir.iter().map(|v| v.abs()).collect::<Vec<_>>());
        trace.peaks.push((*pos, peak));
        let time = peak as f64 / SR_F;
        for (r, &f) in reference.iter_mut().zip(f) {
            *r *= C::from_polar(1.0, 0.0 + (2.0 * PI * f) * time);
        }
        let power: Vec<_> = reference
            .iter()
            .map(|&h| {
                let m = magnitude(h);
                m * m
            })
            .collect();
        let sm = smooth::smooth_power(f, &power, Fraction::PerBin(&frac)).expect("design grid");
        let mp = minimum_phase(&sm.iter().map(|v| v.sqrt() + 1e-12).collect::<Vec<_>>())
            .expect("magnitude");
        for i in 0..BINS {
            let ex = divide(reference[i], mp[i] + C::new(1e-30, 0.0));
            u[i] += divide(C::new(w, 0.0) * ex, C::new(magnitude(ex) + 1e-30, 0.0));
        }
    }
    for u in &mut u {
        *u = divide(*u, C::new(total, 0.0));
    }
    Ok(RcPreparation::Fir(FirPreparation {
        peq,
        p,
        m,
        u,
        nb,
        g,
        lo_fir,
        design_db,
        notes,
    }))
}

fn finish_correction(
    input: &RcInput<'_>,
    phase_frac: f64,
    prepared: &RcPreparation,
    fitter: &mut Fitter<'_>,
    trace: &mut Trace,
    cancelled: &dyn Fn() -> bool,
) -> Result<RcOutput, EngineError> {
    check(cancelled)?;
    let RcPreparation::Fir(prepared) = prepared else {
        let RcPreparation::Complete(out) = prepared else {
            unreachable!()
        };
        return Ok(out.clone());
    };
    let f = design_freqs();
    let (s, c) = (input.setup, input.config);
    let mut peq = prepared.peq.clone();
    let mut p = prepared.p.clone();
    let mut m = prepared.m.clone();
    let u = &prepared.u;
    let nb = &prepared.nb;
    let g = &prepared.g;
    let lo_fir = prepared.lo_fir;
    let design_db = prepared.design_db.clone();
    let notes = prepared.notes.clone();
    let phase =
        unwrap(&smooth::smooth_complex(f, u, Fraction::Scalar(phase_frac)).expect("design grid"));
    let tau = smooth(&gradient(&phase), phase_frac);
    let hf: Vec<_> = f
        .iter()
        .zip(&tau)
        .filter(|(f, _)| **f >= 3000.0 && **f <= 6000.0)
        .map(|(_, t)| *t)
        .collect();
    let reference = decay::median(&hf);
    let bands = band_set(s.rc_lo, s.rc_hi, &s.rc_bands);
    let excess: Vec<_> = f
        .iter()
        .zip(&tau)
        .map(|(&f, &tau)| {
            let fade = phase_fade(f, s, &bands);
            (tau - reference) * fade
        })
        .collect();
    let band: Vec<_> = f.iter().map(|&f| phase_error_band(f, lo_fir, s)).collect();
    let realizer = Realizer {
        excess: &excess,
        band: &band,
        lift: s.lift_ms / 1000.0,
        taps: c.taps,
    };
    let (mut h, mut response, rms0, max0) = realizer.realize(0.0, &m);
    trace.ladder.push((0.0, rms0, max0));
    let mut strength = 0.0;
    if s.lift_ms <= 0.0 {
        (h, response) = truncate(&m, c.taps);
    } else {
        for fraction in [0.9, 0.7, 0.5, 0.3] {
            check(cancelled)?;
            let (hh, ff, rms, max) = realizer.realize(fraction, &m);
            trace.ladder.push((fraction, rms, max));
            if rms <= 0.5_f64.max(rms0 + 0.3) && max <= 2_f64.max(max0 + 1.0) {
                h = hh;
                response = ff;
                strength = fraction;
                break;
            }
        }
    }
    if s.peq_bands > 0 {
        let resid: Vec<_> = (0..BINS)
            .map(|i| {
                let r = design_db[i] - db(response[i]);
                if nb[i] { r.min(0.0) } else { r }
            })
            .collect();
        let cg: Vec<_> = g.iter().map(|&g| interp(g, f, &resid)).collect();
        peq = fit_request(
            input,
            g,
            &cg,
            s.peq_bands,
            20.0,
            500.0,
            Some(s.rc_hi),
            fitter,
        );
        check(cancelled)?;
        let nb_float: Vec<_> = nb.iter().map(|b| if *b { 1.0 } else { 0.0 }).collect();
        for (hz, _, gain) in &mut peq {
            if interp(*hz, f, &nb_float) > 0.5 {
                *gain = gain.min(0.0);
            }
        }
        p = response::peq_response(f, &peq);
        let rc: Vec<_> = f
            .iter()
            .map(|&f| f >= 20.0 && in_bands(f, &bands))
            .collect();
        let design_error = |fx: &[C]| {
            let power: Vec<_> = multiply(fx, &p)
                .iter()
                .map(|&h| {
                    let m = magnitude(h);
                    m * m
                })
                .collect();
            let error: Vec<_> = smooth(&power, 24.0)
                .iter()
                .enumerate()
                .filter(|(i, _)| rc[*i])
                .map(|(i, v)| 10.0 * (v + 1e-30).log10() - design_db[i])
                .collect();
            let med = decay::median(&error);
            rms(&error.iter().map(|e| e - med).collect::<Vec<_>>())
        };
        for _ in 0..2 {
            check(cancelled)?;
            let residual: Vec<_> = (0..BINS)
                .map(|i| {
                    if rc[i] {
                        design_db[i] - db(response[i] * p[i])
                    } else {
                        0.0
                    }
                })
                .collect();
            let broad = smooth(&residual, 3.0);
            let mag: Vec<_> = (0..BINS)
                .map(|i| {
                    let fade = if rc[i] {
                        ((f[i] - lo_fir / 2.0) / (lo_fir / 2.0)).clamp(0.0, 1.0)
                    } else {
                        0.0
                    };
                    magnitude(m[i]) * 10_f64.powf(broad[i] * fade / 20.0)
                })
                .collect();
            let m2 = minimum_phase(&mag).expect("magnitude");
            let (h2, f2) = if s.lift_ms > 0.0 {
                let (h, f, _, _) = realizer.realize(strength, &m2);
                (h, f)
            } else {
                truncate(&m2, c.taps)
            };
            if design_error(&f2) >= design_error(&response) {
                break;
            }
            h = h2;
            response = f2;
            m = m2;
        }
    }
    Ok(RcOutput {
        fir: h,
        response: multiply(&response, &p),
        peq,
        phase_strength: strength,
        phase_band: Some(strength),
        design_db,
        notes,
    })
}

pub(crate) fn phase_score(input: &RcInput<'_>, out: &RcOutput) -> f64 {
    let bands = band_set(input.setup.rc_lo, input.setup.rc_hi, &input.setup.rc_bands);
    let g = band_grid(&bands, 20.0, 1200);
    let gd = response::gd_curve(
        &multiply(main_response(input), &out.response),
        design_freqs(),
        &g,
        Fraction::Scalar(24.0),
    )
    .expect("GD grid");
    let med = decay::median(&gd);
    gd.iter().map(|v| v - med).fold(f64::NEG_INFINITY, f64::max)
}

/// Compare phase candidates while sharing their phase-independent preparation.
/// The supplied fitter must return the same result for the same request: an
/// identical preparation is evaluated once, including the PEQ-only path.
/// Sequential replay clients must account for the shared preparation call.
/// Returns `Err(Cancelled)` once `cancelled()` reads true; the PEQ fitter is
/// expected to watch the same token.
pub fn channel(
    input: &RcInput<'_>,
    fitter: &mut Fitter<'_>,
    cancelled: &dyn Fn() -> bool,
) -> Result<ChannelRc, EngineError> {
    let bands = band_set(input.setup.rc_lo, input.setup.rc_hi, &input.setup.rc_bands);
    channel_with_cuts(
        input,
        fitter,
        &mut |h, lo, hi| {
            if input.setup.rc_bands.is_empty() {
                fit::fit_gd_cuts_until(design_freqs(), h, lo, hi, 6, 3.0, cancelled)
            } else {
                fit::fit_gd_cuts_bands_until(design_freqs(), h, &bands, 6, 3.0, cancelled)
            }
        },
        cancelled,
    )
}
type GdFitter<'a> = dyn FnMut(&[C], f64, f64) -> Vec<Peq> + 'a;

pub(crate) fn channel_with_cuts(
    input: &RcInput<'_>,
    fitter: &mut Fitter<'_>,
    cuts: &mut GdFitter<'_>,
    cancelled: &dyn Fn() -> bool,
) -> Result<ChannelRc, EngineError> {
    let requested = input.config.phase_frac;
    let prepared = prepare_correction(input, fitter, &mut Trace::default(), cancelled)?;
    let mut out = finish_correction(
        input,
        requested,
        &prepared,
        fitter,
        &mut Trace::default(),
        cancelled,
    )?;
    let mut used = None;
    if requested > 6.0 {
        let other = finish_correction(
            input,
            6.0,
            &prepared,
            fitter,
            &mut Trace::default(),
            cancelled,
        )?;
        let mut notes = out.notes.clone();
        notes.extend(other.notes.clone());
        if phase_score(input, &other) < phase_score(input, &out) {
            out = other;
            used = Some(6.0);
        } else {
            used = Some(requested);
        }
        out.notes = notes;
    }
    let gd_cuts = if input.tuning.gd_cuts {
        let cuts = cuts(
            &multiply(main_response(input), &out.response),
            input.setup.rc_lo.max(20.0),
            input.setup.rc_hi,
        );
        // A cancelled GD fit returns meaningless cuts; discard them here.
        check(cancelled)?;
        cuts
    } else {
        Vec::new()
    };
    if input.config.low_protect
        && let Some(hz) = crate::low_protect::rolloff(input.comb, input.main_position, None)
    {
        let mut padded = out.fir.clone();
        padded.resize(N_FFT, 0.0);
        let fir = fft::rfft(&padded);
        let fixed = multiply(&fir, &response::peq_response(design_freqs(), &gd_cuts));
        let limited = crate::low_protect::limit(&mut out.peq, &fixed, hz);
        out.response = multiply(&fir, &response::peq_response(design_freqs(), &out.peq));
        out.notes.push(Note::LowProtect {
            role: input.role.to_owned(),
            rolloff_hz: hz,
            hp_hz: None,
            limited,
        });
    }
    let response = if gd_cuts.is_empty() {
        out.response.clone()
    } else {
        out.notes.push(Note::GdCuts {
            role: input.role.to_owned(),
            cuts: gd_cuts.clone(),
        });
        multiply(
            &out.response,
            &response::peq_response(design_freqs(), &gd_cuts),
        )
    };
    Ok(ChannelRc {
        out,
        phase_frac_used: used,
        gd_cuts,
        response,
    })
}

pub fn combine(
    main: &[C],
    align_delay_ms: f64,
    gain_db: f64,
    xo_hz: Option<f64>,
    sub_sum: Option<&[C]>,
    bus_polarity: i8,
    bus_delay_ms: f64,
) -> Vec<C> {
    assert_eq!(main.len(), BINS);
    let delay = measurement::delay_response(align_delay_ms);
    let gain = C::new(10_f64.powf(gain_db / 20.0), 0.0);
    let mut out: Vec<_> = main.iter().zip(delay).map(|(&h, d)| h * d * gain).collect();
    if let (Some(xo), Some(sub)) = (xo_hz.filter(|v| *v != 0.0), sub_sum) {
        assert_eq!(sub.len(), BINS);
        let hp = response::lr4(design_freqs(), xo, "HP");
        let lp = response::lr4(design_freqs(), xo, "LP");
        let delay = measurement::delay_response(bus_delay_ms);
        for i in 0..BINS {
            out[i] =
                out[i] * hp[i] + sub[i] * C::new(f64::from(bus_polarity), 0.0) * delay[i] * lp[i];
        }
    }
    out
}