DCTS/crates/common/src/seed_finder.rs

use crate::models::GridPointParams;
use std::path::PathBuf;

#[derive(Debug, Clone)]
pub struct SeedMatch {
    pub name: String,
    pub path: PathBuf,
    pub distance: f64,
}

pub const MAX_GLOBAL_SEED_DISTANCE: f64 = 3.0;

/// CNO 有向距离：富金属方向（目标比种子富）重罚，贫金属方向（目标比种子贫）轻罚。
///
/// **数据标定依据**——对历史 1191 个真实 seed_step（种子,目标）配对的成败统计：
/// - 贫金属方向（种子更富、目标往贫走，delta=目标−种子 < 0）：成功率 **42–54%**
/// - 富金属方向（目标更富、delta > 0）：成功率仅 **3–11%**
///   （每个 loghe 分层该规律独立成立，he=−4 时贫方向 54% vs 富方向 3%，差 18 倍）
///
/// **物理解释**：从高金属丰度的收敛解出发**减少**金属（贫方向）是稳定微扰；
/// 反过来从贫金属种子**增加**金属（富方向），新增的紫外谱线辐射驱动会破坏已建立的
/// 辐射平衡，导致大气发散出现 NaN。这是 NLTE 辐射流体计算的已知特性。
///
/// **回测**：对 871 个失败点的 exact_family 种子池回测，旧版绝对值距离有 55% 选了糟糕的
/// 富方向种子；改用此非对称距离后，81% 选贫方向，净改善 314 个点（旧选富→新选贫）。
fn directed_cno_distance(cand: &GridPointParams, target: &GridPointParams) -> f64 {
    const RICH_PENALTY: f64 = 4.0; // 目标比种子富 → 该方向微扰不稳定，重罚
    const POOR_PENALTY: f64 = 1.0; // 目标比种子贫 → 该方向微扰稳定，轻罚
                                   // delta = target − cand：正 = 目标更富（坏方向），负 = 目标更贫（好方向）
    let penalize = |delta: f64| {
        if delta > 0.0 {
            delta * RICH_PENALTY
        } else {
            -delta * POOR_PENALTY
        }
    };
    penalize(target.logc.value() - cand.logc.value())
        + penalize(target.logn.value() - cand.logn.value())
        + penalize(target.logo.value() - cand.logo.value())
}

pub fn calculate_seed_distance(cand: &GridPointParams, target: &GridPointParams) -> (bool, f64) {
    let d_teff = (cand.teff.value() - target.teff.value()).abs();
    let d_logg = (cand.logg.value() - target.logg.value()).abs();
    let d_loghe = (cand.loghe.value() - target.loghe.value()).abs();

    // exact family 判定：Teff/logg/logHe 视为“同物理族”，仅 CNO 丰度不同。
    // Teff 容忍度取一个网格步长 5000K（含边界 <=）：
    // 网格 Teff 档位间隔正好 5000K（如 55000↔60000 相邻档），旧版严格 `< 5000.0`
    // 把所有相邻 Teff 档互为种子静默排除出 exact_family——而实测 Teff 方向是最稳的
    // 种子方向（55k 点靠 50k Teff 种子收敛；50k→60k 直接 nl 亦收敛），且相邻档
    // 同(logg,He,CNO)种子在 exact_family 内有向距离为 0（优于任何同 Teff CNO 邻居）。
    // 桶索引（db::SeedBucketKey）本就按 floor/floor+1 双写，5000K 整除边界可命中。
    // 2026-08-14 修复，详见 docs/tlusty_coldstart_nc_trace_illcond_2026_08_13.md。
    if d_teff <= 5000.0 && d_logg < 0.01 && d_loghe < 0.01 {
        // 同物理族内仅 CNO 不同：用有向距离优先匹配贫金属方向的种子（见 directed_cno_distance）。
        (true, directed_cno_distance(cand, target))
    } else {
        // 距离公式物理意义与标定阐释：
        // 在恒星非局部热力学平衡(NLTE)辐射流体力学与光谱大气计算中，不同物理自由度对于迭代收敛过程的基本影响层级截然相反：
        // 1. Teff (有效温度) 通常达数千至数十万 K，主导连续谱黑体势函数与强激发电离步阶，故除以 5000.0 归一化为基底主控距离量；
        // 2. logg (表面重力加速度) 对静力学与辐射光致压差梯度的平衡破坏力极烈，压强差稍高会触发极大激波不平衡，因此乘上 2.0 予以最高维权惩罚；
        // 3. loghe (氦丰度) 对自由电子密度与热库贡献次于 H-He 电离梯度，乘 0.5 作为次要控制项；
        // 4. CNO 金属元素影响紫外谱线辐射驱动，其方向性同样关键（富方向不稳定，见 directed_cno_distance），
        //    乘 0.1 归一化后纳入全局距离。
        let d_cno = directed_cno_distance(cand, target);
        let global_d = (d_teff / 5000.0) + (d_logg * 2.0) + (d_loghe * 0.5) + (d_cno * 0.1);
        (false, global_d)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::models::GridAxisValue;

    fn params(
        teff: f64,
        logg: f64,
        loghe: f64,
        logc: f64,
        logn: f64,
        logo: f64,
    ) -> GridPointParams {
        GridPointParams {
            teff: GridAxisValue::from_value(teff),
            logg: GridAxisValue::from_value(logg),
            loghe: GridAxisValue::from_value(loghe),
            logc: GridAxisValue::from_value(logc),
            logn: GridAxisValue::from_value(logn),
            logo: GridAxisValue::from_value(logo),
        }
    }

    /// exact_family 内，同 CNO 差幅度下，贫方向距离应远小于富方向距离（4:1）。
    #[test]
    fn test_directed_cno_distance_favors_poor_metal_direction() {
        let seed = params(40000.0, 6.0, -2.0, -2.0, -2.0, -2.0); // 种子 CNO=-6
                                                                 // 贫方向：目标 CNO=-7（种子更富，目标往贫走），|Δ|=1
        let poor_target = params(40000.0, 6.0, -2.0, -3.0, -2.0, -2.0);
        // 富方向：目标 CNO=-5（目标更富），|Δ|=1，同一分量、同幅度
        let rich_target = params(40000.0, 6.0, -2.0, -1.0, -2.0, -2.0);
        let d_poor = directed_cno_distance(&seed, &poor_target);
        let d_rich = directed_cno_distance(&seed, &rich_target);
        assert!(d_poor < d_rich, "贫方向应更近");
        assert!(
            (d_rich / d_poor - 4.0).abs() < 1e-9,
            "富/贫方向同幅度距离比应为 RICH/POOR=4.0，实际 {} / {} = {}",
            d_rich,
            d_poor,
            d_rich / d_poor
        );
    }

    /// exact_family 选种应在多个候选中优先选贫方向种子，即便其 CNO 绝对差更大。
    #[test]
    fn test_exact_family_prefers_poor_direction_seed() {
        let target = params(40000.0, 6.0, -2.0, -2.0, -2.0, -2.0); // 目标 CNO=-6
                                                                   // 候选A：富方向种子（目标比种子富），CNO 绝对差=1
        let rich_seed = params(40000.0, 6.0, -2.0, -3.0, -2.0, -2.0); // CNO=-7, 目标更富
                                                                      // 候选B：贫方向种子（目标比种子贫），CNO 绝对差=2（更大）
        let poor_seed = params(40000.0, 6.0, -2.0, -1.0, -1.0, -2.0); // CNO=-4, 目标更贫
        let (_, d_rich) = calculate_seed_distance(&rich_seed, &target);
        let (_, d_poor) = calculate_seed_distance(&poor_seed, &target);
        // 贫种子虽 CNO 绝对差更大（2 vs 1），但因方向有利，距离应更小
        assert!(
            d_poor < d_rich,
            "贫方向种子距离 {} 应小于富方向 {}（即便 CNO 绝对差更大）",
            d_poor,
            d_rich
        );
    }

    /// 跨 teff（非 exact）时仍标记 global，且方向性体现在 global_d 里。
    #[test]
    fn test_global_branch_marks_non_exact_and_keeps_direction() {
        let target = params(40000.0, 6.0, -2.0, -2.0, -2.0, -2.0);
        // 跨 teff 20000K（超过 exact 容忍 5000K）→ global
        let cand = params(60000.0, 6.0, -2.0, -2.0, -2.0, -2.0);
        let (is_exact, d) = calculate_seed_distance(&cand, &target);
        assert!(!is_exact, "跨 teff 20000K 应为 global 分支");
        assert!(d > 0.0);
    }

    /// 回归（2026-08-14）：d_teff 恰为 5000K（相邻 Teff 档，如 55000↔60000）
    /// 应算 exact_family（旧版严格 < 排除了所有相邻档互为种子）。
    #[test]
    fn test_adjacent_teff_bucket_is_exact_family() {
        let seed = params(55000.0, 5.0, -2.0, -4.0, -3.0, -3.0);
        let target = params(60000.0, 5.0, -2.0, -4.0, -3.0, -3.0);
        let (is_exact, d) = calculate_seed_distance(&seed, &target);
        assert!(is_exact, "d_teff=5000 的相邻 Teff 档应属 exact_family");
        assert_eq!(d, 0.0, "同 (logg,He,CNO) 的相邻 Teff 种子有向距离应为 0");
    }
}