autoeq 0.4.24

//! AutoEQ - A library for audio equalization and filter optimization
//! Loss functions and types for AutoEQ optimizer
//!
//! Copyright (C) 2025-2026 Pierre Aubert pierre(at)spinorama(dot)org
//!
//! This program is free software: you can redistribute it and/or modify
//! it under the terms of the GNU General Public License as published by
//! the Free Software Foundation, either version 3 of the License, or
//! (at your option) any later version.
//!
//! This program is distributed in the hope that it will be useful,
//! but WITHOUT ANY WARRANTY; without even the implied warranty of
//! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//! GNU General Public License for more details.
//!
//! You should have received a copy of the GNU General Public License
//! along with this program.  If not, see <https://www.gnu.org/licenses/>.

use crate::Curve;
use crate::cea2034 as score;
use crate::error::{AutoeqError, Result};
use crate::read;
use clap::ValueEnum;
use ndarray::{Array1, Zip};
use num_complex::Complex64;
use std::collections::HashMap;
use std::f64::consts::PI;

pub mod bass_boost;
pub mod enhanced_weights;
pub mod phase_aware;

/// The type of loss function to use during optimization
#[derive(Debug, Clone, Copy, PartialEq, Eq, ValueEnum)]
pub enum LossType {
    /// Flat loss function (minimize deviation from target curve)
    SpeakerFlat,
    /// Flat loss with asymmetric weighting (peaks penalized 2x more than dips)
    /// Use this for room correction where nulls cannot be fixed with EQ
    SpeakerFlatAsymmetric,
    /// Harman/Olive preference-score objective.
    ///
    /// The helper functions return a preference score, and the optimizer wrapper
    /// converts that score into a minimization objective.
    SpeakerScore,
    /// Flat loss function (minimize deviation from target curve)
    HeadphoneFlat,
    /// Harman headphone preference-score objective.
    ///
    /// The helper functions return a preference score, and the optimizer wrapper
    /// converts that score into a minimization objective.
    HeadphoneScore,
    /// Multi-driver crossover optimization (flatten combined response)
    DriversFlat,
    /// Multi-subwoofer optimization (flatten summed response)
    MultiSubFlat,
    /// EPA (Evaluation, Potency, Activity) perceptual loss.
    /// Combines spectral flatness with sharpness, roughness, and loudness
    /// balance penalties derived from Zwicker psychoacoustic metrics.
    Epa,
}

/// Data required for computing speaker score-based loss
#[derive(Debug, Clone)]
pub struct SpeakerLossData {
    /// On-axis SPL measurements
    pub on: Array1<f64>,
    /// Listening window SPL measurements
    pub lw: Array1<f64>,
    /// Sound power SPL measurements
    pub sp: Array1<f64>,
    /// Predicted in-room SPL measurements
    pub pir: Array1<f64>,
}

impl SpeakerLossData {
    /// Create a new SpeakerLossData instance.
    ///
    /// # Arguments
    /// * `spin` - Map of CEA2034 curves by name ("On Axis", "Listening Window", "Sound Power", "Estimated In-Room Response")
    ///
    /// # Errors
    ///
    /// Returns `AutoeqError::MissingCea2034Curve` if any required curve is missing.
    /// Returns `AutoeqError::CurveLengthMismatch` if curves have different lengths.
    pub fn try_new(spin: &HashMap<String, Curve>) -> Result<Self> {
        let on = spin
            .get("On Axis")
            .ok_or_else(|| AutoeqError::MissingCea2034Curve {
                curve_name: "On Axis".to_string(),
            })?
            .spl
            .clone();
        let lw = spin
            .get("Listening Window")
            .ok_or_else(|| AutoeqError::MissingCea2034Curve {
                curve_name: "Listening Window".to_string(),
            })?
            .spl
            .clone();
        let sp = spin
            .get("Sound Power")
            .ok_or_else(|| AutoeqError::MissingCea2034Curve {
                curve_name: "Sound Power".to_string(),
            })?
            .spl
            .clone();
        let pir = spin
            .get("Estimated In-Room Response")
            .ok_or_else(|| AutoeqError::MissingCea2034Curve {
                curve_name: "Estimated In-Room Response".to_string(),
            })?
            .spl
            .clone();

        // Verify all arrays have the same length
        if on.len() != lw.len() || on.len() != sp.len() || on.len() != pir.len() {
            return Err(AutoeqError::CurveLengthMismatch {
                on_len: on.len(),
                lw_len: lw.len(),
                sp_len: sp.len(),
                pir_len: pir.len(),
            });
        }

        Ok(Self { on, lw, sp, pir })
    }
}

/// Data required for computing headphone loss
#[derive(Debug, Clone)]
pub struct HeadphoneLossData {
    /// Enable smoothing (regularization) of the inverted target curve
    pub smooth: bool,
    /// Smoothing level as 1/N octave (N in [1..24])
    pub smooth_n: usize,
}

impl HeadphoneLossData {
    /// Create a new HeadphoneLossData instance
    ///
    /// # Arguments
    /// * `smooth` - Enable smoothing
    /// * `smooth_n` - Smoothing level as 1/N octave
    pub fn new(smooth: bool, smooth_n: usize) -> Self {
        Self { smooth, smooth_n }
    }
}

/// Crossover filter type for multi-driver optimization
#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
pub enum CrossoverType {
    /// 2nd order Butterworth (12 dB/octave)
    Butterworth2,
    /// 2nd order Linkwitz-Riley (12 dB/octave)
    LinkwitzRiley2,
    /// 4th order Linkwitz-Riley (24 dB/octave) — most common
    #[serde(alias = "LR24")]
    LinkwitzRiley4,
    /// 8th order Linkwitz-Riley (48 dB/octave)
    #[serde(alias = "LR48")]
    LinkwitzRiley8,
    /// No crossover filter (for multi-sub optimization)
    None,
}

impl Default for CrossoverType {
    fn default() -> Self {
        CrossoverType::LinkwitzRiley4
    }
}

impl CrossoverType {
    /// Convert to a short plugin-compatible string.
    pub fn to_plugin_string(&self) -> &'static str {
        match self {
            CrossoverType::Butterworth2 => "Butterworth12",
            CrossoverType::LinkwitzRiley2 => "LR12",
            CrossoverType::LinkwitzRiley4 => "LR24",
            CrossoverType::LinkwitzRiley8 => "LR48",
            CrossoverType::None => "None",
        }
    }

    /// Human-readable display name.
    pub fn display_name(&self) -> &'static str {
        match self {
            CrossoverType::Butterworth2 => "2nd order Butterworth",
            CrossoverType::LinkwitzRiley2 => "2nd order Linkwitz-Riley",
            CrossoverType::LinkwitzRiley4 => "4th order Linkwitz-Riley",
            CrossoverType::LinkwitzRiley8 => "8th order Linkwitz-Riley",
            CrossoverType::None => "No Crossover (Multi-Sub)",
        }
    }
}

impl std::str::FromStr for CrossoverType {
    type Err = String;

    fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
        match s.to_lowercase().as_str() {
            "butterworth2" | "bw2" | "butterworth12" | "bw12" => Ok(CrossoverType::Butterworth2),
            "lr2" | "lr12" | "linkwitzriley2" | "linkwitzriley12" => Ok(CrossoverType::LinkwitzRiley2),
            "lr4" | "lr24" | "linkwitzriley4" | "linkwitzriley24" => Ok(CrossoverType::LinkwitzRiley4),
            "lr8" | "lr48" | "linkwitzriley8" | "linkwitzriley48" => Ok(CrossoverType::LinkwitzRiley8),
            "none" => Ok(CrossoverType::None),
            _ => Err(format!("Unknown crossover type: {}", s)),
        }
    }
}

impl std::fmt::Display for CrossoverType {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.write_str(self.to_plugin_string())
    }
}

/// Measurement data for a single driver
#[derive(Debug, Clone)]
pub struct DriverMeasurement {
    /// Frequency points in Hz
    pub freq: Array1<f64>,
    /// SPL measurements in dB
    pub spl: Array1<f64>,
    /// Phase measurements in degrees (optional for now)
    pub phase: Option<Array1<f64>>,
}

impl DriverMeasurement {
    /// Create a new DriverMeasurement
    pub fn new(freq: Array1<f64>, spl: Array1<f64>, phase: Option<Array1<f64>>) -> Self {
        assert_eq!(freq.len(), spl.len(), "freq and spl must have same length");
        if let Some(ref p) = phase {
            assert_eq!(freq.len(), p.len(), "freq and phase must have same length");
        }
        Self { freq, spl, phase }
    }

    /// Get the frequency range covered by this driver
    pub fn freq_range(&self) -> (f64, f64) {
        let min_freq = self.freq.iter().copied().fold(f64::INFINITY, f64::min);
        let max_freq = self.freq.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        (min_freq, max_freq)
    }

    /// Get the mean frequency (geometric mean)
    pub fn mean_freq(&self) -> f64 {
        let (min_freq, max_freq) = self.freq_range();
        (min_freq * max_freq).sqrt()
    }
}

/// Data required for multi-driver crossover optimization
#[derive(Debug, Clone)]
pub struct DriversLossData {
    /// Measurements for each driver (sorted by frequency range, lowest first)
    pub drivers: Vec<DriverMeasurement>,
    /// Crossover type to use between driver pairs
    pub crossover_type: CrossoverType,
    /// Common frequency grid for evaluation
    pub freq_grid: Array1<f64>,
}

impl DriversLossData {
    /// Create a new DriversLossData instance
    ///
    /// # Arguments
    /// * `drivers` - Vector of driver measurements (will be sorted by frequency)
    /// * `crossover_type` - Type of crossover filter to use
    pub fn new(mut drivers: Vec<DriverMeasurement>, crossover_type: CrossoverType) -> Self {
        assert!(
            drivers.len() >= 2 && drivers.len() <= 4,
            "Must have 2-4 drivers, got {}",
            drivers.len()
        );

        // Sort drivers by their mean frequency (woofer -> midrange -> tweeter)
        drivers.sort_by(|a, b| {
            a.mean_freq()
                .partial_cmp(&b.mean_freq())
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        // Create a common frequency grid spanning all drivers
        // Use logarithmic spacing from lowest to highest frequency
        let min_freq = drivers
            .iter()
            .map(|d| d.freq_range().0)
            .fold(f64::INFINITY, f64::min);
        let max_freq = drivers
            .iter()
            .map(|d| d.freq_range().1)
            .fold(f64::NEG_INFINITY, f64::max);

        // Create log-spaced frequency grid (10 points per octave)
        let freq_grid = crate::read::create_log_frequency_grid(
            10 * 10, // 10 octaves * 10 points per octave
            min_freq.max(20.0),
            max_freq.min(20000.0),
        );

        Self {
            drivers,
            crossover_type,
            freq_grid,
        }
    }
}

/// Compute the flat (current) loss within a specified frequency range
///
/// # Arguments
/// * `freqs` - Frequency points in Hz
/// * `error` - Error values at each frequency point
/// * `min_freq` - Minimum frequency in Hz (inclusive)
/// * `max_freq` - Maximum frequency in Hz (inclusive)
///
/// # Returns
/// * Weighted error value computed only for frequencies within [min_freq, max_freq]
pub fn flat_loss(freqs: &Array1<f64>, error: &Array1<f64>, min_freq: f64, max_freq: f64) -> f64 {
    weighted_mse(freqs, error, min_freq, max_freq)
}

/// Compute the speaker preference score for a candidate PEQ response.
/// `peq_response` must be computed for the candidate parameters.
pub fn speaker_score_loss(
    score_data: &SpeakerLossData,
    freq: &Array1<f64>,
    peq_response: &Array1<f64>,
) -> f64 {
    // Compute 1/2-octave intervals on the fly using the provided frequency grid
    let intervals = score::octave_intervals(2, freq);
    let metrics = if peq_response.iter().all(|v| v.abs() < 1e-12) {
        // Exact score when no PEQ is applied
        score::score(
            freq,
            &intervals,
            &score_data.on,
            &score_data.lw,
            &score_data.sp,
            &score_data.pir,
        )
    } else {
        score::score_peq_approx(
            freq,
            &intervals,
            &score_data.lw,
            &score_data.sp,
            &score_data.pir,
            &score_data.on,
            peq_response,
        )
    };

    metrics.pref_score
}

/// Compute a mixed loss based on flatness on lw and pir
pub fn mixed_loss(
    score_data: &SpeakerLossData,
    freq: &Array1<f64>,
    peq_response: &Array1<f64>,
) -> f64 {
    let lw2 = &score_data.lw + peq_response;
    let pir2 = &score_data.pir + peq_response;
    // Compute slopes in dB per octave over 100 Hz .. 10 kHz
    let lw2_slope = regression_slope_per_octave_in_range(freq, &lw2, 100.0, 10000.0);
    let pir_og_slope = regression_slope_per_octave_in_range(freq, &score_data.pir, 100.0, 10000.0);
    let pir2_slope = regression_slope_per_octave_in_range(freq, &pir2, 100.0, 10000.0);
    if let (Some(lw2eq), Some(pir2og), Some(pir2eq)) = (lw2_slope, pir_og_slope, pir2_slope) {
        // some nlopt algorithms stop for negative values; keep result positive-ish
        (0.5 + lw2eq).powi(2) + (pir2og - pir2eq).powi(2)
    } else {
        f64::INFINITY
    }
}

/// Default bass/treble split frequency in Hz.
///
/// This is a heuristic split point - bass frequencies are weighted more heavily
/// because room acoustics are typically more problematic in the low frequencies.
/// The value of 3000 Hz is based on practical experience but is somewhat arbitrary.
const DEFAULT_BASS_TREBLE_SPLIT_HZ: f64 = 3000.0;

/// Compute weighted mean squared error with frequency-dependent weighting within a frequency range
///
/// # Arguments
/// * `freqs` - Frequency points in Hz
/// * `error` - Error values at each frequency point
/// * `min_freq` - Minimum frequency in Hz (inclusive)
/// * `max_freq` - Maximum frequency in Hz (inclusive)
///
/// # Returns
/// * Weighted error value computed only for frequencies within [min_freq, max_freq]
///
/// # Details
/// Filters frequencies to the specified range, then computes RMS error separately
/// for frequencies below and above the bass/treble split (default 3000 Hz),
/// with higher weight given to the lower frequency band.
/// If the frequency range excludes all data points, returns 0.0.
fn weighted_mse(freqs: &Array1<f64>, error: &Array1<f64>, min_freq: f64, max_freq: f64) -> f64 {
    weighted_mse_with_split(
        freqs,
        error,
        min_freq,
        max_freq,
        DEFAULT_BASS_TREBLE_SPLIT_HZ,
    )
}

/// Compute weighted mean squared error with configurable bass/treble split frequency
///
/// # Arguments
/// * `freqs` - Frequency points in Hz
/// * `error` - Error values at each frequency point
/// * `min_freq` - Minimum frequency in Hz (inclusive)
/// * `max_freq` - Maximum frequency in Hz (inclusive)
/// * `bass_treble_split_hz` - Frequency in Hz that divides bass and treble bands
///
/// # Returns
/// * Weighted error value computed only for frequencies within [min_freq, max_freq]
fn weighted_mse_with_split(
    freqs: &Array1<f64>,
    error: &Array1<f64>,
    min_freq: f64,
    max_freq: f64,
    bass_treble_split_hz: f64,
) -> f64 {
    // Create masks for frequency bands using ndarray's vectorized operations
    let _in_range = freqs.mapv(|f| f >= min_freq && f <= max_freq);
    let bass_band = freqs.mapv(|f| f < bass_treble_split_hz && f >= min_freq && f <= max_freq);
    let treble_band = freqs.mapv(|f| f >= bass_treble_split_hz && f >= min_freq && f <= max_freq);

    // Count points in each band
    let n1: usize = bass_band.iter().filter(|&&b| b).count();
    let n2: usize = treble_band.iter().filter(|&&b| b).count();

    if n1 == 0 && n2 == 0 {
        return f64::INFINITY;
    }

    // Compute squared errors only for valid points
    let squared_errors = error.mapv(|e| e * e);

    let ss1: f64 = Zip::from(&bass_band)
        .and(&squared_errors)
        .fold(0.0, |acc, &mask, &err| if mask { acc + err } else { acc });

    let ss2: f64 = Zip::from(&treble_band)
        .and(&squared_errors)
        .fold(0.0, |acc, &mask, &err| if mask { acc + err } else { acc });

    let err1 = if n1 > 0 {
        (ss1 / n1 as f64).sqrt()
    } else {
        0.0
    };
    let err2 = if n2 > 0 {
        (ss2 / n2 as f64).sqrt()
    } else {
        0.0
    };
    match (n1 > 0, n2 > 0) {
        (true, true) => err1 + err2 / 3.0,
        (true, false) => err1,
        (false, true) => err2,
        (false, false) => f64::INFINITY,
    }
}

/// Configuration for asymmetric loss weighting
///
/// Peaks (positive deviations from target) should be penalized more than dips
/// because:
/// - We can hear peaks clearly and fix them with EQ cuts
/// - Dips/nulls are often caused by acoustic cancellation and cannot be fixed
///   with EQ boosts (boosting into a null just wastes amplifier power)
#[derive(Debug, Clone, Copy)]
pub struct AsymmetricLossConfig {
    /// Weight for positive errors (peaks above transition_freq) - default: 2.0
    pub peak_weight: f64,
    /// Weight for negative errors (dips above transition_freq) - default: 1.0
    pub dip_weight: f64,
    /// Weight for bass peaks (below transition_freq) - default: 5.0
    pub bass_peak_weight: f64,
    /// Weight for bass dips (below transition_freq) - default: 0.2
    pub bass_dip_weight: f64,
    /// Transition frequency between bass and mid/treble weighting - default: 300.0 Hz
    pub transition_freq: f64,
}

impl Default for AsymmetricLossConfig {
    fn default() -> Self {
        Self {
            peak_weight: 2.0,
            dip_weight: 1.0,
            bass_peak_weight: 5.0,
            bass_dip_weight: 0.2,
            transition_freq: 300.0,
        }
    }
}

/// Compute asymmetric weighted mean squared error
///
/// This loss function penalizes peaks (positive errors) more heavily than dips
/// (negative errors), which aligns with psychoacoustic principles:
/// - Peaks are audible and can be corrected with EQ cuts
/// - Dips/nulls cannot be effectively corrected with EQ boosts
///
/// # Arguments
/// * `freqs` - Frequency points in Hz
/// * `error` - Error values at each frequency point (positive = peak, negative = dip)
/// * `min_freq` - Minimum frequency in Hz (inclusive)
/// * `max_freq` - Maximum frequency in Hz (inclusive)
/// * `config` - Asymmetric weighting configuration
///
/// # Returns
/// * Asymmetrically weighted error value
///
/// # Details
/// For each error value:
/// - If error > 0 (peak/overshoot): weighted by `peak_weight`
/// - If error < 0 (dip/undershoot): weighted by `dip_weight`
///
/// Default config uses peak_weight=2.0, dip_weight=1.0 (peaks penalized 2x more)
pub fn weighted_mse_asymmetric(
    freqs: &Array1<f64>,
    error: &Array1<f64>,
    min_freq: f64,
    max_freq: f64,
    config: &AsymmetricLossConfig,
) -> f64 {
    weighted_mse_asymmetric_with_split(
        freqs,
        error,
        min_freq,
        max_freq,
        config,
        DEFAULT_BASS_TREBLE_SPLIT_HZ,
    )
}

/// Compute asymmetric weighted mean squared error with configurable bass/treble split
pub fn weighted_mse_asymmetric_with_split(
    freqs: &Array1<f64>,
    error: &Array1<f64>,
    min_freq: f64,
    max_freq: f64,
    config: &AsymmetricLossConfig,
    bass_treble_split_hz: f64,
) -> f64 {
    // Create masks for frequency bands
    let bass_band = freqs.mapv(|f| f < bass_treble_split_hz && f >= min_freq && f <= max_freq);
    let treble_band = freqs.mapv(|f| f >= bass_treble_split_hz && f >= min_freq && f <= max_freq);

    // Count points in each band
    let n1: usize = bass_band.iter().filter(|&&b| b).count();
    let n2: usize = treble_band.iter().filter(|&&b| b).count();

    if n1 == 0 && n2 == 0 {
        return f64::INFINITY;
    }

    // Compute asymmetrically weighted squared errors with frequency-dependent weights.
    // Below transition_freq: bass_peak_weight / bass_dip_weight
    // Above transition_freq: peak_weight / dip_weight
    // Smooth sigmoid crossfade over ~1 octave centered at transition_freq.
    let log_transition = config.transition_freq.ln();
    // Steepness: ~90% transition within +/-0.5 octave => k = 2*ln(9)/ln(2) ~ 6.34
    let sigmoid_k = 2.0 * 9.0_f64.ln() / 2.0_f64.ln();

    let weighted_squared_errors = Zip::from(freqs).and(error).map_collect(|&f, &e| {
        // sigmoid blend: 0 = full bass weights, 1 = full mid/treble weights
        let blend = 1.0 / (1.0 + (-(f.ln() - log_transition) * sigmoid_k).exp());
        let peak_w =
            config.bass_peak_weight + blend * (config.peak_weight - config.bass_peak_weight);
        let dip_w =
            config.bass_dip_weight + blend * (config.dip_weight - config.bass_dip_weight);
        let weight = if e > 0.0 { peak_w } else { dip_w };
        weight * e * e
    });

    let ss1: f64 = Zip::from(&bass_band)
        .and(&weighted_squared_errors)
        .fold(0.0, |acc, &mask, &err| if mask { acc + err } else { acc });

    let ss2: f64 = Zip::from(&treble_band)
        .and(&weighted_squared_errors)
        .fold(0.0, |acc, &mask, &err| if mask { acc + err } else { acc });

    let err1 = if n1 > 0 {
        (ss1 / n1 as f64).sqrt()
    } else {
        0.0
    };
    let err2 = if n2 > 0 {
        (ss2 / n2 as f64).sqrt()
    } else {
        0.0
    };
    match (n1 > 0, n2 > 0) {
        (true, true) => err1 + err2 / 3.0,
        (true, false) => err1,
        (false, true) => err2,
        (false, false) => f64::INFINITY,
    }
}

/// Compute flat loss with asymmetric weighting (peaks penalized more than dips)
///
/// This is the asymmetric version of `flat_loss()`. Use this when you want the
/// optimizer to prioritize reducing peaks over filling dips.
///
/// # Arguments
/// * `freqs` - Frequency points in Hz
/// * `error` - Error values (positive = peak, negative = dip)
/// * `min_freq` - Minimum frequency in Hz
/// * `max_freq` - Maximum frequency in Hz
///
/// # Returns
/// * Loss value with peaks penalized 2x more than dips
pub fn flat_loss_asymmetric(
    freqs: &Array1<f64>,
    error: &Array1<f64>,
    min_freq: f64,
    max_freq: f64,
) -> f64 {
    weighted_mse_asymmetric(
        freqs,
        error,
        min_freq,
        max_freq,
        &AsymmetricLossConfig::default(),
    )
}

/// Compute the slope (per octave) using linear regression of y against log2(f).
///
/// - `freq`: frequency array in Hz
/// - `y`: corresponding values (e.g., SPL in dB)
/// - Range is defined in Hz as [fmin, fmax]; only f > 0 are considered
/// - Returns `Some(slope_db_per_octave)` or `None` if insufficient data
pub fn regression_slope_per_octave_in_range(
    freq: &Array1<f64>,
    y: &Array1<f64>,
    fmin: f64,
    fmax: f64,
) -> Option<f64> {
    assert_eq!(freq.len(), y.len(), "freq and y must have same length");
    if fmax <= fmin {
        return None;
    }

    let mut n: usize = 0;
    let mut sum_x = 0.0;
    let mut sum_y = 0.0;
    let mut sum_xy = 0.0;
    let mut sum_x2 = 0.0;

    for i in 0..freq.len() {
        let f = freq[i];
        if f > 0.0 && f >= fmin && f <= fmax {
            let xi = f.log2();
            let yi = y[i];
            n += 1;
            sum_x += xi;
            sum_y += yi;
            sum_xy += xi * yi;
            sum_x2 += xi * xi;
        }
    }

    if n < 2 {
        return None;
    }
    let n_f = n as f64;
    let cov_xy = sum_xy - (sum_x * sum_y) / n_f;
    let var_x = sum_x2 - (sum_x * sum_x) / n_f;
    if var_x.abs() < 1e-10 {
        return None;
    }
    Some(cov_xy / var_x)
}

/// Convenience wrapper for slope per octave on a `Curve`.
pub fn curve_slope_per_octave_in_range(curve: &crate::Curve, fmin: f64, fmax: f64) -> Option<f64> {
    regression_slope_per_octave_in_range(&curve.freq, &curve.spl, fmin, fmax)
}

/// Helper to compute complex response of a biquad filter
fn biquad_complex_response(biquad: &crate::iir::Biquad, f: f64) -> Complex64 {
    let (a1, a2, b0, b1, b2) = biquad.constants();
    let omega = 2.0 * PI * f / biquad.srate;
    // z^-1 = e^(-j*omega) = cos(-omega) + j*sin(-omega)
    let z_inv = Complex64::from_polar(1.0, -omega);
    let z_inv2 = z_inv * z_inv;

    let num = b0 + b1 * z_inv + b2 * z_inv2;
    let den = 1.0 + a1 * z_inv + a2 * z_inv2;

    num / den
}

/// Interpolate and prepare driver curves on a common frequency grid
fn prepare_driver_curves(data: &DriversLossData, crossover_freqs: &[f64]) -> Vec<Curve> {
    let n_drivers = data.drivers.len();
    let mut driver_curves = Vec::new();
    for (i, driver) in data.drivers.iter().enumerate() {
        let (passband_low, passband_high) = if let CrossoverType::None = data.crossover_type {
            (20.0, 20000.0)
        } else {
            (
                if i == 0 { 20.0 } else { crossover_freqs[i - 1] },
                if i == n_drivers - 1 {
                    20000.0
                } else {
                    crossover_freqs[i]
                },
            )
        };

        let interpolated = crate::read::normalize_and_interpolate_response_with_range(
            &data.freq_grid,
            &Curve {
                freq: driver.freq.clone(),
                spl: driver.spl.clone(),
                phase: driver.phase.clone(),
            },
            passband_low,
            passband_high,
        );
        driver_curves.push(interpolated);
    }
    driver_curves
}

/// Build crossover filters for a single driver
fn build_crossover_filters_for_driver(
    driver_index: usize,
    n_drivers: usize,
    crossover_type: CrossoverType,
    crossover_freqs: &[f64],
    sample_rate: f64,
) -> Vec<(f64, crate::iir::Biquad)> {
    use crate::iir::{
        peq_butterworth_highpass, peq_butterworth_lowpass, peq_linkwitzriley_highpass,
        peq_linkwitzriley_lowpass,
    };

    let mut filters = Vec::new();
    if let CrossoverType::None = crossover_type {
        return filters;
    }

    if driver_index > 0 {
        let xover_freq = crossover_freqs[driver_index - 1];
        let hp_peq = match crossover_type {
            CrossoverType::Butterworth2 => peq_butterworth_highpass(2, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley2 => peq_linkwitzriley_highpass(2, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley4 => peq_linkwitzriley_highpass(4, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley8 => peq_linkwitzriley_highpass(8, xover_freq, sample_rate),
            CrossoverType::None => vec![],
        };
        filters.extend(hp_peq);
    }
    if driver_index < n_drivers - 1 {
        let xover_freq = crossover_freqs[driver_index];
        let lp_peq = match crossover_type {
            CrossoverType::Butterworth2 => peq_butterworth_lowpass(2, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley2 => peq_linkwitzriley_lowpass(2, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley4 => peq_linkwitzriley_lowpass(4, xover_freq, sample_rate),
            CrossoverType::LinkwitzRiley8 => peq_linkwitzriley_lowpass(8, xover_freq, sample_rate),
            CrossoverType::None => vec![],
        };
        filters.extend(lp_peq);
    }
    filters
}

/// Compute the complex response of a single driver at all frequency points
fn compute_single_driver_complex(
    freq_grid: &Array1<f64>,
    curve: &Curve,
    gain: f64,
    delay_s: f64,
    filters: &[(f64, crate::iir::Biquad)],
) -> Array1<Complex64> {
    let mag_factor = 10.0_f64.powf(gain / 20.0);
    let mut result = Array1::<Complex64>::zeros(freq_grid.len());

    for j in 0..freq_grid.len() {
        let f = freq_grid[j];
        let spl = curve.spl[j];

        let z_driver = if let Some(phase) = &curve.phase {
            let phi = phase[j].to_radians();
            let m = 10.0_f64.powf(spl / 20.0);
            Complex64::from_polar(m, phi)
        } else {
            let m = 10.0_f64.powf(spl / 20.0);
            Complex64::new(m, 0.0)
        };

        let phi_delay = -2.0 * PI * f * delay_s;
        let z_delay = Complex64::from_polar(1.0, phi_delay);

        let mut z_filters = Complex64::new(1.0, 0.0);
        for (_, biquad) in filters {
            z_filters *= biquad_complex_response(biquad, f);
        }

        result[j] = z_driver * mag_factor * z_filters * z_delay;
    }

    result
}

/// Validate driver arguments (shared by combined and per-driver functions)
fn validate_driver_args(
    data: &DriversLossData,
    gains: &[f64],
    crossover_freqs: &[f64],
    delays: Option<&[f64]>,
) {
    let n_drivers = data.drivers.len();
    assert_eq!(gains.len(), n_drivers);
    if !matches!(data.crossover_type, CrossoverType::None) {
        assert_eq!(crossover_freqs.len(), n_drivers - 1);
    }
    if let Some(d) = delays {
        assert_eq!(d.len(), n_drivers);
    }
}

/// Compute the combined frequency response of multiple drivers with crossovers, gains, and delays
///
/// # Arguments
/// * `data` - DriversLossData containing driver measurements and crossover type
/// * `gains` - Gain in dB for each driver
/// * `crossover_freqs` - Crossover frequencies between successive driver pairs
/// * `delays` - Optional delay in ms for each driver
/// * `sample_rate` - Sample rate for filter design
///
/// # Returns
/// * Combined frequency response in dB on the common frequency grid
pub fn compute_drivers_combined_response(
    data: &DriversLossData,
    gains: &[f64],
    crossover_freqs: &[f64],
    delays: Option<&[f64]>,
    sample_rate: f64,
) -> Array1<f64> {
    validate_driver_args(data, gains, crossover_freqs, delays);

    let n_drivers = data.drivers.len();
    let driver_curves = prepare_driver_curves(data, crossover_freqs);

    // Sum complex responses
    let mut combined_complex = Array1::<Complex64>::zeros(data.freq_grid.len());

    for i in 0..n_drivers {
        let delay_s = delays.map(|d| d[i]).unwrap_or(0.0) / 1000.0;
        let filters = build_crossover_filters_for_driver(
            i,
            n_drivers,
            data.crossover_type,
            crossover_freqs,
            sample_rate,
        );
        let driver_complex = compute_single_driver_complex(
            &data.freq_grid,
            &driver_curves[i],
            gains[i],
            delay_s,
            &filters,
        );
        combined_complex += &driver_complex;
    }

    // Convert back to dB SPL
    combined_complex.mapv(|z| 20.0 * z.norm().max(1e-12).log10())
}

/// Compute per-driver frequency responses with crossovers, gains, and delays applied
///
/// Same logic as `compute_drivers_combined_response()` but returns individual
/// driver responses instead of summing them into a single combined response.
///
/// # Arguments
/// * `data` - DriversLossData containing driver measurements and crossover type
/// * `gains` - Gain in dB for each driver
/// * `crossover_freqs` - Crossover frequencies between successive driver pairs
/// * `delays` - Optional delay in ms for each driver
/// * `sample_rate` - Sample rate for filter design
///
/// # Returns
/// * Vec of per-driver frequency responses in dB on the common frequency grid
pub fn compute_per_driver_responses(
    data: &DriversLossData,
    gains: &[f64],
    crossover_freqs: &[f64],
    delays: Option<&[f64]>,
    sample_rate: f64,
) -> Vec<Array1<f64>> {
    validate_driver_args(data, gains, crossover_freqs, delays);

    let n_drivers = data.drivers.len();
    let driver_curves = prepare_driver_curves(data, crossover_freqs);

    let mut results = Vec::with_capacity(n_drivers);
    for i in 0..n_drivers {
        let delay_s = delays.map(|d| d[i]).unwrap_or(0.0) / 1000.0;
        let filters = build_crossover_filters_for_driver(
            i,
            n_drivers,
            data.crossover_type,
            crossover_freqs,
            sample_rate,
        );
        let driver_complex = compute_single_driver_complex(
            &data.freq_grid,
            &driver_curves[i],
            gains[i],
            delay_s,
            &filters,
        );
        results.push(driver_complex.mapv(|z| 20.0 * z.norm().max(1e-12).log10()));
    }

    results
}

/// Compute the loss for multi-driver crossover optimization
///
/// # Arguments
/// * `data` - DriversLossData containing driver measurements
/// * `gains` - Gain in dB for each driver
/// * `crossover_freqs` - Crossover frequencies between successive driver pairs
/// * `sample_rate` - Sample rate for filter design
/// * `min_freq` - Minimum frequency for loss evaluation
/// * `max_freq` - Maximum frequency for loss evaluation
///
/// # Returns
/// * Loss value (lower is better)
pub fn drivers_flat_loss(
    data: &DriversLossData,
    gains: &[f64],
    crossover_freqs: &[f64],
    delays: Option<&[f64]>,
    sample_rate: f64,
    min_freq: f64,
    max_freq: f64,
) -> f64 {
    // Compute combined response
    let combined_response =
        compute_drivers_combined_response(data, gains, crossover_freqs, delays, sample_rate);

    // Normalize the response (subtract the mean in the evaluation range)
    let mut sum = 0.0;
    let mut count = 0;
    for i in 0..data.freq_grid.len() {
        let freq = data.freq_grid[i];
        if freq >= min_freq && freq <= max_freq {
            sum += combined_response[i];
            count += 1;
        }
    }
    let mean = if count > 0 { sum / count as f64 } else { 0.0 };
    let normalized = &combined_response - mean;

    // Compute flatness loss (RMS deviation from zero)
    flat_loss(&data.freq_grid, &normalized, min_freq, max_freq)
}

/// Compute the loss for multi-subwoofer optimization (flat summed response)
///
/// # Arguments
/// * `data` - DriversLossData
/// * `gains` - Gain in dB for each sub
/// * `delays` - Delay in ms for each sub
/// * `sample_rate` - Sample rate
/// * `min_freq` - Min freq for evaluation
/// * `max_freq` - Max freq for evaluation
///
/// # Returns
/// * Loss value
pub fn multisub_flat_loss(
    data: &DriversLossData,
    gains: &[f64],
    delays: &[f64],
    sample_rate: f64,
    min_freq: f64,
    max_freq: f64,
) -> f64 {
    // Pass empty crossover freqs (ignored because CrossoverType::None)
    let crossover_freqs = vec![];
    let combined_response =
        compute_drivers_combined_response(data, gains, &crossover_freqs, Some(delays), sample_rate);

    // Normalize the response (subtract the mean in the evaluation range)
    let mut sum = 0.0;
    let mut count = 0;
    for i in 0..data.freq_grid.len() {
        let freq = data.freq_grid[i];
        if freq >= min_freq && freq <= max_freq {
            sum += combined_response[i];
            count += 1;
        }
    }
    let mean = if count > 0 { sum / count as f64 } else { 0.0 };
    let normalized = &combined_response - mean;

    // Compute flatness loss (RMS deviation from zero)
    flat_loss(&data.freq_grid, &normalized, min_freq, max_freq)
}

/// Calculate the standard deviation (SD) of the deviation error over the specified frequency range.
///
/// This function filters the input curve to include only frequencies within the specified range,
/// then calculates the standard deviation of the deviation values.
///
/// # Arguments
/// * `freq` - Frequency array in Hz
/// * `deviation` - Deviation values from Harman target curve in dB
/// * `fmin` - Minimum frequency in Hz (typically 50 Hz)
/// * `fmax` - Maximum frequency in Hz (typically 10000 Hz)
///
/// # Returns
/// * Standard deviation of the deviation in the specified frequency range
///
/// # Notes
/// Used as part of the Olive et al. headphone preference prediction model.
pub fn calculate_standard_deviation_in_range(
    freq: &Array1<f64>,
    deviation: &Array1<f64>,
    fmin: f64,
    fmax: f64,
) -> f64 {
    assert_eq!(
        freq.len(),
        deviation.len(),
        "freq and deviation must have same length"
    );

    let mut values = Vec::new();

    // Collect deviation values in the specified frequency range
    for i in 0..freq.len() {
        let f = freq[i];
        if f >= fmin && f <= fmax {
            values.push(deviation[i]);
        }
    }

    if values.is_empty() {
        return 0.0;
    }

    // Calculate mean
    let mean = values.iter().sum::<f64>() / values.len() as f64;

    // Calculate variance
    let variance = values.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / values.len() as f64;

    // Return standard deviation
    variance.sqrt()
}

/// Calculate the absolute slope (AS) of the deviation using logarithmic regression over the specified frequency range.
///
/// This function performs linear regression of deviation against log2(frequency) to determine
/// the slope, then returns the absolute value.
///
/// # Arguments
/// * `freq` - Frequency array in Hz
/// * `deviation` - Deviation values from Harman target curve in dB
/// * `fmin` - Minimum frequency in Hz (typically 50 Hz)
/// * `fmax` - Maximum frequency in Hz (typically 10000 Hz)
///
/// # Returns
/// * Absolute value of the slope in dB per octave
///
/// # Notes
/// Used as part of the Olive et al. headphone preference prediction model.
fn calculate_absolute_slope_in_range(
    freq: &Array1<f64>,
    deviation: &Array1<f64>,
    fmin: f64,
    fmax: f64,
) -> f64 {
    match regression_slope_per_octave_in_range(freq, deviation, fmin, fmax) {
        Some(slope) => slope.abs(),
        None => 0.0,
    }
}

/// Compute headphone preference score based on frequency response deviations
///
/// This implements the headphone preference prediction model from:
/// Olive, S. E., Welti, T., & McMullin, E. (2013). "A Statistical Model that
/// Predicts Listeners' Preference Ratings of Around-Ear and On-Ear Headphones"
///
/// The model predicts preference using the equation:
/// **Predicted Preference Rating = 114.49 - (12.62 × SD) - (15.52 × AS)**
///
/// Where:
/// - **SD** = Standard deviation of the deviation error over 50 Hz to 10 kHz
/// - **AS** = Absolute value of slope of the deviation over 50 Hz to 10 kHz
///
/// # Arguments
/// * `curve` - Frequency response curve representing deviation from Harman AE/OE target
///
/// # Returns
/// * Predicted preference rating (higher values indicate better preference)
///
/// # Important Note
/// The input curve should represent deviation from the Harman Around-Ear (AE) or
/// On-Ear (OE) target curve, **not** deviation from flat or neutral response.
///
/// The frequency range for calculations is 50 Hz to 10 kHz as specified in the paper.
///
/// # References
/// Olive, S. E., Welti, T., & McMullin, E. (2013). "A Statistical Model that
/// Predicts Listeners' Preference Ratings of Around-Ear and On-Ear Headphones".
/// Presented at the 135th Convention of the Audio Engineering Society.
pub fn headphone_loss(curve: &Curve) -> f64 {
    let freq = &curve.freq;
    let deviation = &curve.spl;

    // Define frequency range for analysis (50 Hz to 10 kHz per paper)
    const FMIN: f64 = 50.0;
    const FMAX: f64 = 10000.0;

    // Calculate SD (Standard Deviation) of the deviation error
    let sd = calculate_standard_deviation_in_range(freq, deviation, FMIN, FMAX);

    // Calculate AS (Absolute Slope) of the deviation
    let as_value = calculate_absolute_slope_in_range(freq, deviation, FMIN, FMAX);

    // Apply the Olive et al. equation (Equation 4 from the paper)
    // Predicted Preference Rating = 114.49 - (12.62 × SD) - (15.52 × AS)

    // Return the preference rating directly.
    // Optimizer wrappers convert this score into a minimization objective.
    114.49 - (12.62 * sd) - (15.52 * as_value)
}

/// Compute headphone preference score with additional target curve
///
/// # Arguments
/// * `data` - Headphone loss data containing smoothing parameters
/// * `response` - Measured frequency response in dB
/// * `target` - Target frequency response in dB
///
/// # Returns
/// * Predicted headphone preference score where higher is better
pub fn headphone_loss_with_target(
    data: &HeadphoneLossData,
    response: &Curve,
    target: &Curve,
) -> f64 {
    // freqs on which we normalize every curve: 12 points per octave between 20 and 20kHz
    let freqs = read::create_log_frequency_grid(10 * 12, 20.0, 20000.0);

    let input_curve = read::normalize_and_interpolate_response(&freqs, response);
    let target_curve = read::normalize_and_interpolate_response(&freqs, target);

    // normalized and potentially smooth
    let deviation = Curve {
        freq: freqs.clone(),
        spl: &target_curve.spl - &input_curve.spl,
        phase: None,
    };
    let smooth_deviation = if data.smooth {
        read::smooth_one_over_n_octave(&deviation, data.smooth_n)
    } else {
        deviation.clone()
    };

    headphone_loss(&smooth_deviation)
}

#[cfg(test)]
mod tests {
    use super::*;
    use ndarray::Array1;
    use ndarray::array;
    use std::collections::HashMap;

    #[test]
    fn score_loss_matches_score_when_peq_zero() {
        // Simple synthetic data
        let freq = Array1::from(vec![100.0, 1000.0]);
        let on = Array1::from(vec![80.0_f64, 85.0_f64]);
        let lw = Array1::from(vec![81.0_f64, 84.0_f64]);
        let sp = Array1::from(vec![78.0_f64, 82.0_f64]);
        let pir = Array1::from(vec![80.5_f64, 84.0_f64]);

        // Build spin map expected by constructor
        let mut spin: HashMap<String, Curve> = HashMap::new();
        spin.insert(
            "On Axis".to_string(),
            Curve {
                freq: freq.clone(),
                spl: on.clone(),
                phase: None,
            },
        );
        spin.insert(
            "Listening Window".to_string(),
            Curve {
                freq: freq.clone(),
                spl: lw.clone(),
                phase: None,
            },
        );
        spin.insert(
            "Sound Power".to_string(),
            Curve {
                freq: freq.clone(),
                spl: sp.clone(),
                phase: None,
            },
        );
        spin.insert(
            "Estimated In-Room Response".to_string(),
            Curve {
                freq: freq.clone(),
                spl: pir.clone(),
                phase: None,
            },
        );

        let sd = SpeakerLossData::try_new(&spin).expect("test spin data should be valid");
        let zero = Array1::zeros(freq.len());

        // Expected preference using score() with zero PEQ (i.e., base curves)
        let intervals = super::score::octave_intervals(2, &freq);
        let expected = super::score::score(&freq, &intervals, &on, &lw, &sp, &pir);
        let got = speaker_score_loss(&sd, &freq, &zero);
        if got.is_nan() && expected.pref_score.is_nan() {
            // ok
        } else {
            assert!((got - expected.pref_score).abs() < 1e-12);
        }
    }

    #[test]
    fn regression_slope_per_octave_linear_log_relation_full_range() {
        // y = 3 * log2(f) + 1
        let freq = Array1::from(vec![100.0, 200.0, 400.0, 800.0]);
        let y = freq.mapv(|f: f64| 3.0 * f.log2() + 1.0);
        let slope = regression_slope_per_octave_in_range(&freq, &y, 100.0, 800.0).unwrap();
        assert!((slope - 3.0).abs() < 1e-12);
    }

    #[test]
    fn regression_slope_per_octave_sub_range() {
        // Same log-linear relation, sub-range 200..=800
        let freq = Array1::from(vec![100.0, 200.0, 400.0, 800.0]);
        let y = freq.mapv(|f: f64| -2.5 * f.log2() + 4.0);
        let slope = regression_slope_per_octave_in_range(&freq, &y, 200.0, 800.0).unwrap();
        assert!((slope + 2.5).abs() < 1e-12);
    }

    #[test]
    fn test_calculate_standard_deviation_in_range() {
        // Test SD calculation with known values
        let freq = Array1::from(vec![50.0, 100.0, 1000.0, 5000.0, 10000.0]);
        let deviation = Array1::from(vec![1.0, 2.0, 3.0, 4.0, 5.0]); // All in range

        let sd = calculate_standard_deviation_in_range(&freq, &deviation, 50.0, 10000.0);

        // Manual calculation: mean = (1+2+3+4+5)/5 = 3.0
        // variance = ((1-3)² + (2-3)² + (3-3)² + (4-3)² + (5-3)²)/5 = (4+1+0+1+4)/5 = 2.0
        // sd = sqrt(2.0) ≈ 1.414
        let expected_sd = 2.0_f64.sqrt();
        assert!(
            (sd - expected_sd).abs() < 1e-12,
            "SD calculation incorrect: got {}, expected {}",
            sd,
            expected_sd
        );
    }

    #[test]
    fn test_calculate_standard_deviation_filtered_range() {
        // Test SD calculation with frequency filtering
        let freq = Array1::from(vec![20.0, 100.0, 1000.0, 5000.0, 15000.0]); // Some out of range
        let deviation = Array1::from(vec![10.0, 2.0, 4.0, 6.0, 20.0]); // First and last should be filtered

        let sd = calculate_standard_deviation_in_range(&freq, &deviation, 50.0, 10000.0);

        // Only values at 100Hz, 1kHz, 5kHz should be included: [2.0, 4.0, 6.0]
        // mean = (2+4+6)/3 = 4.0
        // variance = ((2-4)² + (4-4)² + (6-4)²)/3 = (4+0+4)/3 = 8/3
        // sd = sqrt(8/3) ≈ 1.633
        let expected_sd = (8.0_f64 / 3.0_f64).sqrt();
        assert!(
            (sd - expected_sd).abs() < 1e-12,
            "SD calculation with filtering incorrect: got {}, expected {}",
            sd,
            expected_sd
        );
    }

    #[test]
    fn test_calculate_absolute_slope_in_range() {
        // Test AS calculation with linear slope
        let freq = Array1::from(vec![
            50.0, 100.0, 200.0, 400.0, 800.0, 1600.0, 3200.0, 6400.0, 10000.0,
        ]);
        // Create a perfect 2 dB/octave slope: y = 2 * log2(f) + constant
        let deviation = freq.mapv(|f: f64| 2.0 * f.log2());

        let as_value = calculate_absolute_slope_in_range(&freq, &deviation, 50.0, 10000.0);

        // Should return absolute value of 2.0
        assert!(
            (as_value - 2.0).abs() < 1e-12,
            "AS calculation incorrect: got {}, expected 2.0",
            as_value
        );
    }

    #[test]
    fn test_calculate_absolute_slope_negative() {
        // Test AS calculation with negative slope
        let freq = Array1::from(vec![
            50.0, 100.0, 200.0, 400.0, 800.0, 1600.0, 3200.0, 6400.0, 10000.0,
        ]);
        // Create a perfect -3 dB/octave slope
        let deviation = freq.mapv(|f: f64| -3.0 * f.log2());

        let as_value = calculate_absolute_slope_in_range(&freq, &deviation, 50.0, 10000.0);

        // Should return absolute value of -3.0 = 3.0
        assert!(
            (as_value - 3.0).abs() < 1e-12,
            "AS calculation with negative slope incorrect: got {}, expected 3.0",
            as_value
        );
    }

    #[test]
    fn test_headphone_loss_perfect_harman_deviation() {
        // Test with zero deviation from Harman target (perfect response)
        let freq = Array1::from(vec![50.0, 100.0, 1000.0, 5000.0, 10000.0]);
        let deviation = Array1::zeros(5); // Perfect match to Harman target

        let curve = Curve {
            freq: freq.clone(),
            spl: deviation,
            phase: None,
        };
        let score = headphone_loss(&curve);

        // With zero deviation (SD=0, AS=0), predicted preference = 114.49
        // The helper returns the raw preference score.
        let expected_score = 114.49;
        assert!(
            (score - expected_score).abs() < 1e-12,
            "Perfect Harman score incorrect: got {}, expected {}",
            score,
            expected_score
        );
    }

    #[test]
    fn test_headphone_loss_with_deviation() {
        // Test with some deviation from Harman target
        let freq = Array1::from(vec![50.0, 100.0, 1000.0, 5000.0, 10000.0]);
        let deviation = Array1::from(vec![1.0, 1.0, 1.0, 1.0, 1.0]); // 1dB constant deviation

        let curve = Curve {
            freq: freq.clone(),
            spl: deviation,
            phase: None,
        };
        let score = headphone_loss(&curve);

        // SD = 0 (constant deviation), AS = 0 (flat slope)
        // predicted preference = 114.49 - (12.62 * 0) - (15.52 * 0) = 114.49
        // But wait - SD should be 0 for constant values, but AS should also be 0
        let expected_preference = 114.49;
        let expected_score = expected_preference;
        assert!(
            (score - expected_score).abs() < 1e-10,
            "Constant deviation score incorrect: got {}, expected {}",
            score,
            expected_score
        );
    }

    #[test]
    fn test_headphone_loss_with_slope() {
        // Test with a sloped deviation
        let freq = Array1::from(vec![
            50.0, 100.0, 200.0, 400.0, 800.0, 1600.0, 3200.0, 6400.0, 10000.0,
        ]);
        // Create a 1 dB/octave slope in the deviation
        let deviation = freq.mapv(|f: f64| 1.0 * f.log2());

        let curve = Curve {
            freq: freq.clone(),
            spl: deviation,
            phase: None,
        };
        let score = headphone_loss(&curve);

        // AS = 1.0 (absolute slope)
        // SD will be non-zero due to the slope
        // predicted preference = 114.49 - (12.62 * SD) - (15.52 * 1.0)
        // Score should be worse (more negative) than perfect case (-114.49)
        assert!(
            score > 50.0,
            "Sloped deviation should have lower preference: got {}",
            score
        );
    }

    #[test]
    fn test_headphone_loss_with_target() {
        // Test the target curve variant with zero deviation
        let freq = Array1::logspace(10.0, 1.301, 4.301, 100);
        let response = Array1::from_elem(100, 5.0); // Constant 5dB response
        let target = Array1::from_elem(100, 5.0); // Same as response

        let response_curve = Curve {
            freq: freq.clone(),
            spl: response,
            phase: None,
        };
        let target_curve = Curve {
            freq: freq.clone(),
            spl: target,
            phase: None,
        };
        let data = HeadphoneLossData::new(false, 2);
        let score = headphone_loss_with_target(&data, &response_curve, &target_curve);

        // When response matches target, deviation is zero, so should get perfect score
        let expected_perfect_score = 114.49;
        assert!(
            (score - expected_perfect_score).abs() < 1e-10,
            "Perfect target match score incorrect: got {}, expected {}",
            score,
            expected_perfect_score
        );
    }

    #[test]
    fn mixed_loss_finite_with_zero_peq() {
        // Frequency grid
        let freq = Array1::from(vec![
            100.0, 200.0, 400.0, 800.0, 1600.0, 3200.0, 6400.0, 10000.0,
        ]);
        // Zero curves
        let on = Array1::zeros(freq.len());
        let lw = Array1::zeros(freq.len());
        let sp = Array1::zeros(freq.len());
        let pir = Array1::zeros(freq.len());

        // Build spin map
        let mut spin: HashMap<String, Curve> = HashMap::new();
        spin.insert(
            "On Axis".to_string(),
            Curve {
                freq: freq.clone(),
                spl: on,
                phase: None,
            },
        );
        spin.insert(
            "Listening Window".to_string(),
            Curve {
                freq: freq.clone(),
                spl: lw,
                phase: None,
            },
        );
        spin.insert(
            "Sound Power".to_string(),
            Curve {
                freq: freq.clone(),
                spl: sp,
                phase: None,
            },
        );
        spin.insert(
            "Estimated In-Room Response".to_string(),
            Curve {
                freq: freq.clone(),
                spl: pir,
                phase: None,
            },
        );

        let sd = SpeakerLossData::try_new(&spin).expect("test spin data should be valid");
        let peq = Array1::zeros(freq.len());
        let v = mixed_loss(&sd, &freq, &peq);
        assert!(v.is_finite(), "mixed_loss should be finite, got {}", v);
    }

    #[test]
    fn test_weighted_mse_basic() {
        // Two points below 3k, two points above 3k with unit error
        let freqs = array![1000.0, 2000.0, 4000.0, 8000.0];
        let err = array![1.0, 1.0, 1.0, 1.0];
        let v = weighted_mse(&freqs, &err, 100.0, 10000.0); // Full range
        // RMS below = 1, RMS above = 1 -> total = 1 + 1/3 = 1.333...
        assert!((v - (1.0 + 1.0 / 3.0)).abs() < 1e-12, "got {}", v);
    }

    #[test]
    fn test_weighted_mse_empty_upper_segment() {
        // All freqs below 3k -> upper RMS = 0
        let freqs = array![100.0, 200.0, 500.0];
        let err = array![2.0, 2.0, 2.0]; // squares: 4,4,4 -> mean=4 -> rms=2
        let v = weighted_mse(&freqs, &err, 50.0, 10000.0); // Full range
        assert!((v - 2.0).abs() < 1e-12, "got {}", v);
    }

    #[test]
    fn test_weighted_mse_scaling() {
        // Different errors below and above to verify weighting
        let freqs = array![1000.0, 1500.0, 4000.0, 6000.0];
        let err = array![2.0, 2.0, 3.0, 3.0];
        // below RMS = sqrt((4+4)/2)=2, above RMS = sqrt((9+9)/2)=3
        let v = weighted_mse(&freqs, &err, 500.0, 10000.0); // Full range
        let expected = 2.0 + 3.0 / 3.0; // 3.0
        assert!((v - expected).abs() < 1e-12, "got {}", v);
    }

    #[test]
    fn test_weighted_mse_frequency_filtering() {
        // Test that frequency filtering works correctly
        let freqs = array![100.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0];
        let err = array![1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 1.0];

        // Filter to only include 1kHz-4kHz range
        let v = weighted_mse(&freqs, &err, 1000.0, 4000.0);
        // Only frequencies 1000, 2000, 4000 should be included with errors 2, 2, 3
        // All below 3kHz: RMS = sqrt((4+4)/2) = 2
        // 4kHz above 3kHz: RMS = sqrt(9/1) = 3
        let expected = 2.0 + 3.0 / 3.0; // 3.0
        assert!(
            (v - expected).abs() < 1e-12,
            "got {} expected {}",
            v,
            expected
        );
    }

    #[test]
    fn test_weighted_mse_no_frequencies_in_range() {
        // Test edge case where no frequencies are in range
        let freqs = array![100.0, 200.0, 500.0];
        let err = array![2.0, 3.0, 1.0];

        // Filter to range that excludes all frequencies
        let v = weighted_mse(&freqs, &err, 1000.0, 5000.0);
        assert!(
            v.is_infinite(),
            "Should return INFINITY when no frequencies in range"
        );
    }

    #[test]
    fn test_weighted_mse_partial_range() {
        // Test filtering that includes only some frequencies
        let freqs = array![100.0, 1000.0, 2000.0, 4000.0, 8000.0];
        let err = array![1.0, 2.0, 2.0, 3.0, 4.0];

        // Filter to 500-3000 range (should include 1000, 2000)
        let v = weighted_mse(&freqs, &err, 500.0, 3000.0);
        // Only 1000, 2000 with errors 2, 2 - both below 3kHz
        // RMS = sqrt((4+4)/2) = 2.0, no high freq component
        let expected = 2.0;
        assert!(
            (v - expected).abs() < 1e-12,
            "got {} expected {}",
            v,
            expected
        );
    }

    #[test]
    fn test_flat_loss_frequency_filtering() {
        // Test that flat_loss correctly delegates to weighted_mse with frequency bounds
        let freqs = array![100.0, 1000.0, 2000.0, 4000.0, 8000.0];
        let err = array![1.0, 2.0, 2.0, 3.0, 4.0];

        let v1 = flat_loss(&freqs, &err, 1000.0, 4000.0);
        let v2 = weighted_mse(&freqs, &err, 1000.0, 4000.0);

        assert_eq!(v1, v2, "flat_loss should equal weighted_mse");
    }

    #[test]
    fn test_frequency_filtering_boundary_conditions() {
        // Test inclusive boundaries
        let freqs = array![1000.0, 2000.0, 3000.0];
        let err = array![1.0, 1.0, 1.0];

        // Include only exact boundary frequencies
        let v = weighted_mse(&freqs, &err, 1000.0, 3000.0);
        // All three frequencies should be included
        // 1000, 2000 are below 3kHz threshold (n1=2, ss1=2)
        // 3000 is >= 3kHz threshold (n2=1, ss2=1)
        // err1 = sqrt(2/2) = 1.0, err2 = sqrt(1/1) = 1.0
        // result = 1.0 + 1.0/3.0 = 1.333...
        let expected = 1.0 + 1.0 / 3.0;
        assert!(
            (v - expected).abs() < 1e-12,
            "got {} expected {}",
            v,
            expected
        );

        // Exclude boundary frequencies
        let v2 = weighted_mse(&freqs, &err, 1001.0, 2999.0);
        // Only 2000 Hz should be included (below 3kHz threshold)
        // err1 = sqrt(1/1) = 1.0, err2 = 0
        // result = 1.0 + 0/3.0 = 1.0
        let expected2 = 1.0;
        assert!(
            (v2 - expected2).abs() < 1e-12,
            "got {} expected {}",
            v2,
            expected2
        );
    }

    #[test]
    fn test_headphone_loss_perfect_correction() {
        // Test that zero deviation gives perfect score
        let freq = Array1::logspace(10.0, 1.699, 4.0, 100); // 50Hz to 10kHz
        let zero_deviation = Array1::zeros(100);

        let curve = Curve {
            freq: freq.clone(),
            spl: zero_deviation,
            phase: None,
        };
        let score = headphone_loss(&curve);

        // Perfect correction should give score of 114.49
        let expected_perfect = 114.49;
        assert!(
            (score - expected_perfect).abs() < 1e-10,
            "Perfect correction score incorrect: got {}, expected {}",
            score,
            expected_perfect
        );
    }

    #[test]
    fn test_headphone_loss_sign_independence() {
        // Test that headphone_loss gives same result for +deviation and -deviation
        // (since SD and AS are sign-independent)
        let freq = Array1::logspace(10.0, 1.699, 4.0, 100);

        // Create a deviation with varying values
        let deviation_positive = freq.mapv(|f: f64| 0.5 * f.log2() + 2.0);
        let deviation_negative = -&deviation_positive;

        let curve_pos = Curve {
            freq: freq.clone(),
            spl: deviation_positive,
            phase: None,
        };
        let curve_neg = Curve {
            freq: freq.clone(),
            spl: deviation_negative,
            phase: None,
        };

        let score_pos = headphone_loss(&curve_pos);
        let score_neg = headphone_loss(&curve_neg);

        // Scores should be equal since SD is symmetric and AS uses absolute value
        assert!(
            (score_pos - score_neg).abs() < 1e-10,
            "Sign independence violated: pos={}, neg={}",
            score_pos,
            score_neg
        );
    }

    #[test]
    fn test_headphone_loss_worse_than_perfect() {
        // Test that non-zero deviation gives worse score than zero
        let freq = Array1::logspace(10.0, 1.699, 4.0, 100);
        let zero_deviation = Array1::zeros(100);
        let nonzero_deviation = Array1::from_elem(100, 3.0); // 3dB constant deviation

        let perfect_curve = Curve {
            freq: freq.clone(),
            spl: zero_deviation,
            phase: None,
        };
        let imperfect_curve = Curve {
            freq: freq.clone(),
            spl: nonzero_deviation,
            phase: None,
        };

        let perfect_score = headphone_loss(&perfect_curve);
        let imperfect_score = headphone_loss(&imperfect_curve);

        // Imperfect should score lower (worse) than perfect
        assert!(
            imperfect_score < perfect_score,
            "Imperfect correction should score lower: perfect={}, imperfect={}",
            perfect_score,
            imperfect_score
        );

        // Perfect should be 114.49
        assert!(
            (perfect_score - 114.49).abs() < 1e-10,
            "Perfect score should be 114.49, got {}",
            perfect_score
        );
    }

    #[test]
    fn test_weighted_mse_treble_only() {
        // All frequencies above 3kHz — bass band is empty
        let freqs = Array1::from_vec(vec![4000.0, 8000.0, 16000.0]);
        let error = Array1::from_vec(vec![2.0, 2.0, 2.0]);
        let result = weighted_mse(&freqs, &error, 4000.0, 16000.0);
        // With all treble, should return err2 (not err2/3)
        assert!(
            (result - 2.0).abs() < 1e-10,
            "treble-only loss should be full RMS, got {}",
            result
        );
    }

    #[test]
    fn test_weighted_mse_bass_only() {
        // All frequencies below 3kHz — treble band is empty
        let freqs = Array1::from_vec(vec![100.0, 500.0, 2000.0]);
        let error = Array1::from_vec(vec![3.0, 3.0, 3.0]);
        let result = weighted_mse(&freqs, &error, 100.0, 2000.0);
        // With all bass, should return err1
        assert!(
            (result - 3.0).abs() < 1e-10,
            "bass-only loss should be full RMS, got {}",
            result
        );
    }

    #[test]
    fn test_weighted_mse_both_bands() {
        // Mix of bass and treble
        let freqs = Array1::from_vec(vec![100.0, 1000.0, 5000.0, 10000.0]);
        let error = Array1::from_vec(vec![3.0, 3.0, 6.0, 6.0]);
        let result = weighted_mse(&freqs, &error, 100.0, 10000.0);
        // err1 (bass RMS) = 3.0, err2 (treble RMS) = 6.0
        // result = 3.0 + 6.0/3.0 = 5.0
        assert!(
            (result - 5.0).abs() < 1e-10,
            "two-band loss incorrect, got {}",
            result
        );
    }

    #[test]
    fn test_regression_slope_identical_freqs() {
        // All frequencies identical — var_x should be ~0, return None
        let freq = Array1::from_vec(vec![1000.0, 1000.0, 1000.0]);
        let y = Array1::from_vec(vec![80.0, 85.0, 90.0]);
        let result = regression_slope_per_octave_in_range(&freq, &y, 999.0, 1001.0);
        assert!(
            result.is_none(),
            "identical frequencies should return None, got {:?}",
            result
        );
    }

    #[test]
    fn test_bass_asymmetry_penalizes_peaks_heavily() {
        // A 10dB peak at 80Hz should produce ~5x the penalty of same peak at 2kHz
        // because bass_peak_weight=5.0 vs peak_weight=2.0
        let config = AsymmetricLossConfig::default();

        // Single-frequency test at 80Hz (deep bass)
        let freqs_bass = Array1::from_vec(vec![80.0]);
        let error_bass = Array1::from_vec(vec![10.0]); // 10dB peak

        let loss_bass = weighted_mse_asymmetric_with_split(
            &freqs_bass,
            &error_bass,
            20.0,
            20000.0,
            &config,
            3000.0,
        );

        // Single-frequency test at 2kHz (mid range)
        let freqs_mid = Array1::from_vec(vec![2000.0]);
        let error_mid = Array1::from_vec(vec![10.0]); // 10dB peak

        let loss_mid = weighted_mse_asymmetric_with_split(
            &freqs_mid,
            &error_mid,
            20.0,
            20000.0,
            &config,
            3000.0,
        );

        // 80Hz is well below 300Hz transition, so weight ~ bass_peak_weight = 5.0
        // 2kHz is well above 300Hz transition, so weight ~ peak_weight = 2.0
        // The loss returns sqrt(weight * e^2) = e * sqrt(weight)
        // So ratio = sqrt(5.0) / sqrt(2.0) = sqrt(2.5) ~ 1.58
        let ratio = loss_bass / loss_mid;
        let expected_ratio = (5.0_f64 / 2.0).sqrt();
        assert!(
            (ratio - expected_ratio).abs() < 0.1,
            "bass peak penalty ratio should be ~{:.2}x (sqrt(5/2)), got {:.2}",
            expected_ratio,
            ratio
        );
    }

    #[test]
    fn test_bass_dips_nearly_ignored() {
        // 10dB dip at 80Hz should produce near-zero penalty (bass_dip_weight=0.2)
        let config = AsymmetricLossConfig::default();

        // 10dB dip at 80Hz
        let freqs = Array1::from_vec(vec![80.0]);
        let error_dip = Array1::from_vec(vec![-10.0]);

        let loss_dip = weighted_mse_asymmetric_with_split(
            &freqs,
            &error_dip,
            20.0,
            20000.0,
            &config,
            3000.0,
        );

        // 10dB peak at 80Hz for comparison
        let error_peak = Array1::from_vec(vec![10.0]);
        let loss_peak = weighted_mse_asymmetric_with_split(
            &freqs,
            &error_peak,
            20.0,
            20000.0,
            &config,
            3000.0,
        );

        // dip weight 0.2 vs peak weight 5.0, so dip loss / peak loss ~ 0.2/5.0 = 0.04
        let ratio = loss_dip / loss_peak;
        assert!(
            ratio < 0.25,
            "bass dip penalty should be much smaller than bass peak, ratio={:.4}",
            ratio
        );
    }
}