#!/usr/bin/env python3
"""
红宝石压力标定系统 —— 最优拟合算法 v2.0
=======================================
相比原版洛伦兹拟合算法的改进：
  1. 双峰 pseudo-Voigt 替代纯 Lorentzian（兼顾自然线宽 + 仪器展宽）
  2. 多项式自适应基线（解决非均匀背景问题）
  3. 鲁棒初始值估计（自适应 ROI + 二阶导零交叉定位）
  4. 加权拟合（峰区高权重，远离峰的噪声低权重）
  5. Bootstrap 不确定度估计（比协方差矩阵更稳健）
  6. 多次采集自动平均与统计
  7. AIC/BIC 模型选择（自动判断 Lorentzian vs Voigt 哪个更优）
  8. 支持批量文件处理

依赖: numpy, scipy, matplotlib, pandas
安装: pip install numpy scipy matplotlib pandas
"""

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit, minimize
from scipy.signal import savgol_filter, find_peaks
from scipy.special import wofz   # 用于精确 Voigt
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')

# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第一部分：线型函数库                                                ║
# ╚══════════════════════════════════════════════════════════════════════╝

def lorentzian(x, x0, gamma, A):
    """纯 Lorentzian 线型（不含基线）"""
    return A * (gamma / 2) ** 2 / ((x - x0) ** 2 + (gamma / 2) ** 2)


def gaussian(x, x0, sigma, A):
    """纯 Gaussian 线型（不含基线）"""
    fwhm_g = 2.3548 * sigma
    return A * np.exp(-4 * np.log(2) * (x - x0) ** 2 / fwhm_g ** 2)


def voigt_exact(x, x0, sigma, gamma, A):
    """
    精确 Voigt 线型（基于 Faddeeva 函数）
    sigma : Gaussian 分量的标准差
    gamma : Lorentzian 分量的 FWHM
    """
    z = ((x - x0) + 1j * gamma / 2) / (sigma * np.sqrt(2))
    w = wofz(z)
    profile = np.real(w) / (sigma * np.sqrt(2 * np.pi))
    # 归一化使峰高为 A
    z0 = (1j * gamma / 2) / (sigma * np.sqrt(2))
    peak_val = np.real(wofz(z0)) / (sigma * np.sqrt(2 * np.pi))
    return A * profile / peak_val if peak_val > 0 else profile


def pseudo_voigt(x, x0, fwhm, A, eta):
    """
    伪 Voigt 线型
    eta : Lorentzian 分量比例 (0=纯高斯, 1=纯洛伦兹)
    fwhm : 总半高全宽
    """
    eta = np.clip(eta, 0, 1)
    L = (fwhm / 2) ** 2 / ((x - x0) ** 2 + (fwhm / 2) ** 2)
    G = np.exp(-4 * np.log(2) * (x - x0) ** 2 / fwhm ** 2)
    return A * (eta * L + (1 - eta) * G)


def voigt_fwhm(sigma, gamma):
    """
    Voigt 线型的 FWHM 近似公式 (Thompson et al. 1987)
    sigma : Gaussian 标准差
    gamma : Lorentzian FWHM
    """
    fG = 2.3548 * sigma
    fL = gamma
    return (fG**5 + 2.69*fG**4*fL + 2.43*fG**3*fL**2
            + 4.47*fG**2*fL**3 + 0.07*fG*fL**4 + fL**5) ** 0.2


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第二部分：双峰模型                                                  ║
# ╚══════════════════════════════════════════════════════════════════════╝

def double_voigt(x, x1, x2, sigma1, sigma2, gamma1, gamma2, A1, A2, *bg_coeffs):
    """
    双峰精确 Voigt + 多项式基线
    bg_coeffs: 多项式系数 (c0, c1, c2, ...)，基线 = c0 + c1*(x-x_c) + c2*(x-x_c)^2
    """
    x_c = np.mean(x)
    baseline = np.zeros_like(x)
    for i, c in enumerate(bg_coeffs):
        baseline += c * (x - x_c) ** i

    return (baseline
            + voigt_exact(x, x1, sigma1, gamma1, A1)
            + voigt_exact(x, x2, sigma2, gamma2, A2))


def double_pseudo_voigt(x, x1, x2, fwhm1, fwhm2, A1, A2, eta1, eta2, *bg_coeffs):
    """
    双峰伪 Voigt + 多项式基线（计算更快，适合实时场景）
    """
    x_c = np.mean(x)
    baseline = np.zeros_like(x)
    for i, c in enumerate(bg_coeffs):
        baseline += c * (x - x_c) ** i

    return (baseline
            + pseudo_voigt(x, x1, fwhm1, A1, eta1)
            + pseudo_voigt(x, x2, fwhm2, A2, eta2))


def double_lorentzian(x, x1, x2, g1, g2, A1, A2, *bg_coeffs):
    """双峰 Lorentzian + 多项式基线（兼容原版）"""
    x_c = np.mean(x)
    baseline = np.zeros_like(x)
    for i, c in enumerate(bg_coeffs):
        baseline += c * (x - x_c) ** i

    return (baseline
            + lorentzian(x, x1, g1, A1)
            + lorentzian(x, x2, g2, A2))


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第三部分：压力计算                                                  ║
# ╚══════════════════════════════════════════════════════════════════════╝

# 多种常用 Ruby 压力标度
PRESSURE_SCALES = {
    'mao_hydrostatic': {'A': 1904.0, 'B': 7.665, 'R1_0': 694.22,
                        'ref': 'Mao et al. 1986 (hydrostatic)'},
    'mao_nonhydrostatic': {'A': 1904.0, 'B': 5.0, 'R1_0': 694.22,
                           'ref': 'Mao et al. 1986 (non-hydrostatic)'},
    'dewaele': {'A': 1920.0, 'B': 9.61, 'R1_0': 694.24,
                'ref': 'Dewaele et al. 2004'},
    'dorogokupets': {'A': 1876.0, 'B': 10.71, 'R1_0': 694.24,
                     'ref': 'Dorogokupets & Oganov 2007'},
}


def ruby_pressure(lambda_r1_nm, scale='mao_nonhydrostatic'):
    """
    输入 R1 峰位 (nm)，输出压力 (GPa)
    支持多种标度选择
    """
    s = PRESSURE_SCALES[scale]
    A, B, R1_0 = s['A'], s['B'], s['R1_0']
    delta = lambda_r1_nm - R1_0
    P = (A / B) * ((1 + B * delta / R1_0) ** B - 1)
    return P


def ruby_pressure_gradient(lambda_r1_nm, scale='mao_nonhydrostatic'):
    """dP/dλ 用于误差传播"""
    s = PRESSURE_SCALES[scale]
    A, B, R1_0 = s['A'], s['B'], s['R1_0']
    return A * (1 + B * (lambda_r1_nm - R1_0) / R1_0) ** (B - 1) / R1_0


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第四部分：数据加载与预处理                                          ║
# ╚══════════════════════════════════════════════════════════════════════╝

def load_spectrum(filepath):
    """
    加载光谱 CSV 文件
    自动检测元数据行并跳过
    返回: wavelength, intensity, metadata_dict
    """
    filepath = Path(filepath)
    meta = {}

    with open(filepath, 'r') as f:
        lines = f.readlines()

    # 解析元数据（前几行非数据行）
    data_start = 0
    for i, line in enumerate(lines):
        parts = line.strip().split(',')
        if len(parts) == 2:
            try:
                float(parts[0])
                float(parts[1])
                data_start = i
                break
            except ValueError:
                meta[parts[0].strip()] = parts[1].strip()

    df = pd.read_csv(filepath, skiprows=data_start, header=None,
                     names=['Wavelength', 'Intensity'])
    return df['Wavelength'].values, df['Intensity'].values, meta


def estimate_noise(wavelength, intensity, peak_wl, margin=5):
    """
    估计噪声水平
    在远离峰值的基线区域计算标准差
    """
    # 取光谱两端各 margin nm 作为基线区域
    left_mask = wavelength < (wavelength.min() + margin)
    right_mask = wavelength > (wavelength.max() - margin)
    # 排除峰区
    baseline_mask = (left_mask | right_mask) & (np.abs(wavelength - peak_wl) > 3)

    if baseline_mask.sum() < 10:
        # 回退：用最低 10% 的数据估算
        threshold = np.percentile(intensity, 10)
        baseline_mask = intensity < threshold

    return np.std(intensity[baseline_mask])


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第五部分：鲁棒初始值估计                                            ║
# ╚══════════════════════════════════════════════════════════════════════╝

def find_ruby_peaks(wavelength, intensity, sg_window=15, sg_order=2):
    """
    鲁棒的 R1/R2 峰位初始估计
    策略：
      1. S-G 平滑
      2. 自适应 ROI 检测
      3. 二阶导数零交叉辅助定位
      4. 多重回退机制
    """
    # S-G 平滑
    win = min(sg_window, len(intensity) // 4 * 2 + 1)
    if win < 5:
        win = 5
    smoothed = savgol_filter(intensity, window_length=win, polyorder=sg_order)

    # 估计背景和信号
    bg_level = np.percentile(intensity, 10)
    peak_max = np.max(smoothed)
    signal = peak_max - bg_level

    # 全局最大值作为主峰
    main_idx = np.argmax(smoothed)
    main_wl = wavelength[main_idx]

    # 自适应 ROI：主峰 ± 4nm
    roi_half = 4.0
    roi_mask = (wavelength > main_wl - roi_half) & (wavelength < main_wl + roi_half)
    wl_roi = wavelength[roi_mask]
    sm_roi = smoothed[roi_mask]
    int_roi = intensity[roi_mask]

    # 尝试找两个峰
    # 方法1: find_peaks 自适应 prominence
    for prom_frac in [0.02, 0.01, 0.005, 0.05, 0.1]:
        peaks_idx, props = find_peaks(sm_roi, distance=5,
                                       prominence=signal * prom_frac)
        if len(peaks_idx) >= 2:
            break

    if len(peaks_idx) >= 2:
        # 取最高的两个
        top2 = peaks_idx[np.argsort(sm_roi[peaks_idx])[::-1]][:2]
        top2 = np.sort(top2)
        x_R2 = wl_roi[top2[0]]  # 短波长侧 = R2
        x_R1 = wl_roi[top2[1]]  # 长波长侧 = R1
        A_R2 = sm_roi[top2[0]] - bg_level
        A_R1 = sm_roi[top2[1]] - bg_level
    else:
        # 方法2: 二阶导数零交叉
        if len(sm_roi) > 10:
            d2 = np.gradient(np.gradient(sm_roi, wl_roi), wl_roi)
            d2_smooth = savgol_filter(d2, window_length=min(11, len(d2)//2*2+1),
                                      polyorder=2)
            # 二阶导最小值 = 峰位
            neg_peaks, _ = find_peaks(-d2_smooth, distance=5)
            if len(neg_peaks) >= 2:
                top2 = neg_peaks[np.argsort(-sm_roi[neg_peaks])][:2]
                top2 = np.sort(top2)
                x_R2 = wl_roi[top2[0]]
                x_R1 = wl_roi[top2[1]]
                A_R2 = sm_roi[top2[0]] - bg_level
                A_R1 = sm_roi[top2[1]] - bg_level
            else:
                # 方法3: 假设 R1/R2 间距约 1.4nm
                x_R1 = main_wl
                x_R2 = main_wl - 1.4
                A_R1 = signal
                A_R2 = signal * 0.5
        else:
            x_R1 = main_wl
            x_R2 = main_wl - 1.4
            A_R1 = signal
            A_R2 = signal * 0.5

    return {
        'x_R1': x_R1, 'x_R2': x_R2,
        'A_R1': max(A_R1, 1), 'A_R2': max(A_R2, 1),
        'bg': bg_level,
        'roi_mask': roi_mask,
        'wl_roi': wl_roi,
        'int_roi': int_roi,
        'smoothed_roi': sm_roi,
        'signal': signal,
    }


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第六部分：核心拟合引擎                                              ║
# ╚══════════════════════════════════════════════════════════════════════╝

def fit_ruby(wavelength, counts, model='auto', baseline_order=1,
             weighted=True, bootstrap_n=200, scale='mao_nonhydrostatic',
             plot=True, title=''):
    """
    红宝石荧光光谱拟合 —— 最优解

    参数:
    ------
    wavelength : np.array, 波长 (nm)
    counts     : np.array, 强度
    model      : 'voigt', 'pseudo_voigt', 'lorentzian', 'auto'
                 auto 模式下自动通过 AIC 选择最优模型
    baseline_order : 基线多项式阶数 (0=常数, 1=线性, 2=二次)
    weighted   : 是否使用加权拟合
    bootstrap_n: Bootstrap 重采样次数 (0=跳过，用协方差矩阵)
    scale      : 压力标度名称
    plot       : 是否绘图
    title      : 图标题附加信息

    返回:
    ------
    dict: 包含所有拟合结果
    """
    # ── Step 1: 初始值估计 ──
    peaks = find_ruby_peaks(wavelength, counts)
    wl_roi = peaks['wl_roi']
    int_roi = peaks['int_roi']
    x_R1, x_R2 = peaks['x_R1'], peaks['x_R2']
    A_R1, A_R2 = peaks['A_R1'], peaks['A_R2']
    bg = peaks['bg']

    # ── Step 2: 构建权重 ──
    if weighted:
        # 基于信号强度的权重：峰区高权重，基线区低权重
        noise_est = estimate_noise(wavelength, counts, x_R1)
        noise_est = max(noise_est, 1.0)  # 防止除零
        # Poisson-like: sigma ∝ sqrt(counts), 加上读出噪声
        sigma_i = np.sqrt(np.maximum(int_roi, 1) + noise_est ** 2)
        weights = 1.0 / sigma_i
    else:
        weights = np.ones_like(int_roi)

    # ── Step 3: 定义并拟合多个模型 ──
    results = {}
    n_bg = baseline_order + 1
    bg_init = [bg] + [0.0] * baseline_order

    # --- Lorentzian ---
    def _fit_lorentz():
        p0 = [x_R1, x_R2, 0.5, 0.4, A_R1, A_R2] + bg_init
        n_params = len(p0)
        lo = [x_R1-2, x_R2-2, 0.05, 0.05, 0, 0] + [-np.inf]*n_bg
        hi = [x_R1+2, x_R2+2, 5, 5, A_R1*10, A_R2*10] + [np.inf]*n_bg
        lo[6] = 0  # bg >= 0

        popt, pcov = curve_fit(double_lorentzian, wl_roi, int_roi,
                               p0=p0, bounds=(lo, hi), sigma=1/weights,
                               absolute_sigma=False, maxfev=20000)
        fitted = double_lorentzian(wl_roi, *popt)
        return popt, pcov, fitted, n_params, 'Lorentzian'

    # --- Pseudo-Voigt ---
    def _fit_pvoigt():
        p0 = [x_R1, x_R2, 0.5, 0.4, A_R1, A_R2, 0.5, 0.5] + bg_init
        n_params = len(p0)
        lo = [x_R1-2, x_R2-2, 0.05, 0.05, 0, 0, 0, 0] + [-np.inf]*n_bg
        hi = [x_R1+2, x_R2+2, 5, 5, A_R1*10, A_R2*10, 1, 1] + [np.inf]*n_bg
        lo[8] = 0  # bg >= 0

        popt, pcov = curve_fit(double_pseudo_voigt, wl_roi, int_roi,
                               p0=p0, bounds=(lo, hi), sigma=1/weights,
                               absolute_sigma=False, maxfev=20000)
        fitted = double_pseudo_voigt(wl_roi, *popt)
        return popt, pcov, fitted, n_params, 'Pseudo-Voigt'

    # --- Exact Voigt ---
    def _fit_voigt():
        p0 = [x_R1, x_R2, 0.15, 0.12, 0.3, 0.25, A_R1, A_R2] + bg_init
        n_params = len(p0)
        lo = [x_R1-2, x_R2-2, 0.01, 0.01, 0.01, 0.01, 0, 0] + [-np.inf]*n_bg
        hi = [x_R1+2, x_R2+2, 2, 2, 3, 3, A_R1*10, A_R2*10] + [np.inf]*n_bg
        lo[8] = 0  # bg >= 0

        popt, pcov = curve_fit(double_voigt, wl_roi, int_roi,
                               p0=p0, bounds=(lo, hi), sigma=1/weights,
                               absolute_sigma=False, maxfev=20000)
        fitted = double_voigt(wl_roi, *popt)
        return popt, pcov, fitted, n_params, 'Voigt'

    # 尝试所有模型
    fit_funcs = {
        'lorentzian': _fit_lorentz,
        'pseudo_voigt': _fit_pvoigt,
        'voigt': _fit_voigt,
    }

    models_to_try = (['lorentzian', 'pseudo_voigt', 'voigt']
                     if model == 'auto' else [model])

    for m_name in models_to_try:
        try:
            popt, pcov, fitted, n_p, label = fit_funcs[m_name]()
            perr = np.sqrt(np.diag(pcov))

            # AIC / BIC
            n = len(int_roi)
            residuals = int_roi - fitted
            ss_res = np.sum(residuals ** 2)
            ss_tot = np.sum((int_roi - np.mean(int_roi)) ** 2)
            r_squared = 1 - ss_res / ss_tot
            adj_r_squared = 1 - (1 - r_squared) * (n - 1) / (n - n_p - 1)

            log_likelihood = -n/2 * np.log(2*np.pi*ss_res/n) - n/2
            aic = 2 * n_p - 2 * log_likelihood
            bic = n_p * np.log(n) - 2 * log_likelihood

            results[m_name] = {
                'popt': popt, 'pcov': pcov, 'perr': perr,
                'fitted': fitted, 'n_params': n_p, 'label': label,
                'R_squared': r_squared, 'adj_R_squared': adj_r_squared,
                'AIC': aic, 'BIC': bic, 'residuals': residuals,
            }
        except Exception as e:
            results[m_name] = {'error': str(e)}

    # ── Step 4: 模型选择 ──
    valid_models = {k: v for k, v in results.items() if 'error' not in v}
    if not valid_models:
        raise RuntimeError(f"所有模型拟合失败: {results}")

    # 选 AIC 最小的
    best_name = min(valid_models, key=lambda k: valid_models[k]['AIC'])
    best = valid_models[best_name]

    # 提取 R1 峰位（所有模型第一个参数都是 x_R1）
    r1_fit = best['popt'][0]
    r1_err = best['perr'][0]
    r2_fit = best['popt'][1]
    r2_err = best['perr'][1]

    # 提取 FWHM
    if best_name == 'lorentzian':
        fwhm_r1 = best['popt'][2]
        fwhm_r2 = best['popt'][3]
    elif best_name == 'pseudo_voigt':
        fwhm_r1 = best['popt'][2]
        fwhm_r2 = best['popt'][3]
    elif best_name == 'voigt':
        fwhm_r1 = voigt_fwhm(best['popt'][2], best['popt'][4])
        fwhm_r2 = voigt_fwhm(best['popt'][3], best['popt'][5])

    # ── Step 5: Bootstrap 不确定度 ──
    if bootstrap_n > 0:
        r1_bootstrap = []
        r2_bootstrap = []
        rng = np.random.default_rng(42)

        for _ in range(bootstrap_n):
            # 残差重采样
            boot_residuals = rng.choice(best['residuals'], size=len(int_roi),
                                        replace=True)
            boot_data = best['fitted'] + boot_residuals

            try:
                if best_name == 'lorentzian':
                    bp, _ = curve_fit(double_lorentzian, wl_roi, boot_data,
                                      p0=best['popt'], maxfev=5000)
                elif best_name == 'pseudo_voigt':
                    bp, _ = curve_fit(double_pseudo_voigt, wl_roi, boot_data,
                                      p0=best['popt'], maxfev=5000)
                elif best_name == 'voigt':
                    bp, _ = curve_fit(double_voigt, wl_roi, boot_data,
                                      p0=best['popt'], maxfev=5000)
                r1_bootstrap.append(bp[0])
                r2_bootstrap.append(bp[1])
            except Exception:
                pass

        if len(r1_bootstrap) > 10:
            r1_boot_std = np.std(r1_bootstrap)
            r2_boot_std = np.std(r2_bootstrap)
            # 取 Bootstrap 和协方差估计中较大的作为保守估计
            r1_err = max(r1_err, r1_boot_std)
            r2_err = max(r2_err, r2_boot_std)
            bootstrap_success = True
            bootstrap_count = len(r1_bootstrap)
        else:
            bootstrap_success = False
            bootstrap_count = len(r1_bootstrap)
    else:
        bootstrap_success = None
        bootstrap_count = 0

    # ── Step 6: 压力计算 ──
    P = ruby_pressure(r1_fit, scale=scale)
    dPdL = ruby_pressure_gradient(r1_fit, scale=scale)
    P_err = abs(dPdL) * r1_err

    # ── Step 7: SNR ──
    noise_std = estimate_noise(wavelength, counts, r1_fit)
    snr = peaks['signal'] / noise_std if noise_std > 0 else np.nan

    # ── Step 8: 输出 ──
    print("=" * 60)
    print(f"  最优模型 : {best['label']}  (AIC={best['AIC']:.1f})")
    print(f"  R1 峰位  : {r1_fit:.4f} ± {r1_err:.4f} nm")
    print(f"  R2 峰位  : {r2_fit:.4f} ± {r2_err:.4f} nm")
    print(f"  FWHM R1  : {fwhm_r1:.4f} nm")
    print(f"  FWHM R2  : {fwhm_r2:.4f} nm")
    print(f"  压  力   : {P:.3f} ± {P_err:.3f} GPa  ({PRESSURE_SCALES[scale]['ref']})")
    print(f"  R²       : {best['R_squared']:.6f}")
    print(f"  adj-R²   : {best['adj_R_squared']:.6f}")
    print(f"  SNR      : {snr:.1f}" if not np.isnan(snr) else "  SNR      : —")
    if bootstrap_success is not None:
        status = f"成功 ({bootstrap_count}/{bootstrap_n})" if bootstrap_success else "样本不足"
        print(f"  Bootstrap: {status}")
    if len(valid_models) > 1:
        print("  ── 模型对比 ──")
        for mn, mv in sorted(valid_models.items(), key=lambda x: x[1]['AIC']):
            marker = " ◀ BEST" if mn == best_name else ""
            print(f"    {mv['label']:15s}  AIC={mv['AIC']:8.1f}  "
                  f"BIC={mv['BIC']:8.1f}  R²={mv['R_squared']:.6f}{marker}")
    print("=" * 60)

    # ── Step 9: 绘图 ──
    if plot:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        import matplotlib.gridspec as gridspec

        fig = plt.figure(figsize=(12, 8))
        gs = gridspec.GridSpec(2, 2, height_ratios=[3, 1], hspace=0.05, wspace=0.3)

        # 主拟合图
        ax1 = fig.add_subplot(gs[0, 0])
        ax1.plot(wl_roi, int_roi, 'o', ms=1.5, color='#378ADD', alpha=0.5, label='Raw data')

        wl_fine = np.linspace(wl_roi.min(), wl_roi.max(), 1000)
        if best_name == 'lorentzian':
            fit_fine = double_lorentzian(wl_fine, *best['popt'])
            # 分峰
            bg_fine = np.zeros_like(wl_fine)
            x_c = np.mean(wl_fine)
            for i, c in enumerate(best['popt'][6:]):
                bg_fine += c * (wl_fine - x_c) ** i
            p1 = bg_fine + lorentzian(wl_fine, best['popt'][0], best['popt'][2], best['popt'][4])
            p2 = bg_fine + lorentzian(wl_fine, best['popt'][1], best['popt'][3], best['popt'][5])
        elif best_name == 'pseudo_voigt':
            fit_fine = double_pseudo_voigt(wl_fine, *best['popt'])
            bg_fine = np.zeros_like(wl_fine)
            x_c = np.mean(wl_fine)
            for i, c in enumerate(best['popt'][8:]):
                bg_fine += c * (wl_fine - x_c) ** i
            p1 = bg_fine + pseudo_voigt(wl_fine, best['popt'][0], best['popt'][2],
                                         best['popt'][4], best['popt'][6])
            p2 = bg_fine + pseudo_voigt(wl_fine, best['popt'][1], best['popt'][3],
                                         best['popt'][5], best['popt'][7])
        elif best_name == 'voigt':
            fit_fine = double_voigt(wl_fine, *best['popt'])
            bg_fine = np.zeros_like(wl_fine)
            x_c = np.mean(wl_fine)
            for i, c in enumerate(best['popt'][8:]):
                bg_fine += c * (wl_fine - x_c) ** i
            p1 = bg_fine + voigt_exact(wl_fine, best['popt'][0], best['popt'][2],
                                        best['popt'][4], best['popt'][6])
            p2 = bg_fine + voigt_exact(wl_fine, best['popt'][1], best['popt'][3],
                                        best['popt'][5], best['popt'][7])

        ax1.plot(wl_fine, fit_fine, '-', color='#D85A30', lw=2,
                label=f'{best["label"]} fit (R²={best["R_squared"]:.5f})')
        ax1.plot(wl_fine, p1, '--', color='#D85A30', lw=1, alpha=0.6,
                label=f'R1 = {r1_fit:.3f} nm')
        ax1.plot(wl_fine, p2, '--', color='#639922', lw=1, alpha=0.6,
                label=f'R2 = {r2_fit:.3f} nm')
        ax1.plot(wl_fine, bg_fine, ':', color='gray', lw=1, alpha=0.5, label='Baseline')

        ax1.set_ylabel('Intensity (counts)')
        ax1.legend(fontsize=8, loc='upper right')
        ax1.set_title(f'Ruby Fluorescence  |  P = {P:.2f} ± {P_err:.2f} GPa  |  '
                      f'FWHM(R1) = {fwhm_r1:.3f} nm\n{title}',
                      fontsize=10)
        ax1.set_xticklabels([])
        ax1.grid(True, alpha=0.2)

        # 残差图
        ax2 = fig.add_subplot(gs[1, 0], sharex=ax1)
        ax2.plot(wl_roi, best['residuals'], '.', ms=1.5, color='#378ADD', alpha=0.5)
        ax2.axhline(0, color='gray', lw=0.5)
        ax2.set_xlabel('Wavelength (nm)')
        ax2.set_ylabel('Residual')
        ax2.grid(True, alpha=0.2)

        # 模型对比（如果有多个）
        ax3 = fig.add_subplot(gs[0, 1])
        if len(valid_models) > 1:
            names = [v['label'] for v in valid_models.values()]
            aics = [v['AIC'] for v in valid_models.values()]
            colors = ['#D85A30' if n == best['label'] else '#378ADD' for n in names]
            bars = ax3.barh(names, aics, color=colors, alpha=0.8, edgecolor='black', lw=0.5)
            ax3.set_xlabel('AIC (lower is better)')
            ax3.set_title('Model Comparison', fontsize=10)
            for bar, val in zip(bars, aics):
                ax3.text(bar.get_width(), bar.get_y() + bar.get_height()/2,
                        f' {val:.1f}', va='center', fontsize=9)
        else:
            ax3.text(0.5, 0.5, f'Single model:\n{best["label"]}',
                    transform=ax3.transAxes, ha='center', va='center', fontsize=12)
            ax3.set_frame_on(False)
            ax3.set_xticks([])
            ax3.set_yticks([])

        # 残差直方图
        ax4 = fig.add_subplot(gs[1, 1])
        ax4.hist(best['residuals'], bins=30, color='#378ADD', alpha=0.7,
                edgecolor='black', lw=0.5, density=True)
        # 正态分布参考
        x_hist = np.linspace(best['residuals'].min(), best['residuals'].max(), 100)
        res_std = np.std(best['residuals'])
        ax4.plot(x_hist, 1/(res_std*np.sqrt(2*np.pi)) * np.exp(-x_hist**2/(2*res_std**2)),
                'r-', lw=1.5, label='Normal')
        ax4.set_xlabel('Residual')
        ax4.set_ylabel('Density')
        ax4.set_title('Residual Distribution', fontsize=10)
        ax4.legend(fontsize=8)
        ax4.grid(True, alpha=0.2)

        plt.tight_layout()
        fig_path = Path(wavelength.__class__.__module__).parent if hasattr(wavelength, '__module__') else Path('.')
        fig_path = Path('/sessions/zen-blissful-pasteur')
        plt.savefig(fig_path / 'ruby_fit_result.png', dpi=150, bbox_inches='tight')
        plt.close()
        print(f"  图已保存: {fig_path / 'ruby_fit_result.png'}")

    return {
        'model': best_name,
        'model_label': best['label'],
        'R1_nm': r1_fit, 'R1_err': r1_err,
        'R2_nm': r2_fit, 'R2_err': r2_err,
        'FWHM_R1': fwhm_r1, 'FWHM_R2': fwhm_r2,
        'P_GPa': P, 'P_err': P_err,
        'scale': scale,
        'R_squared': best['R_squared'],
        'adj_R_squared': best['adj_R_squared'],
        'AIC': best['AIC'], 'BIC': best['BIC'],
        'SNR': snr,
        'bootstrap_success': bootstrap_success,
        'all_models': {k: {'AIC': v['AIC'], 'BIC': v['BIC'], 'R_squared': v['R_squared']}
                       for k, v in valid_models.items()},
        'popt': best['popt'],
    }


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第七部分：批量处理与多次平均                                        ║
# ╚══════════════════════════════════════════════════════════════════════╝

def batch_fit(file_list, average_first=False, **kwargs):
    """
    批量拟合多个光谱文件

    参数:
    ------
    file_list    : list of str/Path, 光谱文件路径列表
    average_first: True = 先将多个光谱平均再拟合（同条件重复测量推荐）
                   False = 分别拟合每个文件再统计
    **kwargs     : 传给 fit_ruby 的其他参数

    返回:
    ------
    dict: 汇总统计结果
    """
    kwargs.setdefault('plot', False)
    file_list = [Path(f) for f in file_list]

    if average_first and len(file_list) > 1:
        # 加载并对齐波长轴（假设相同仪器配置，波长轴一致）
        all_wl = []
        all_int = []
        for f in file_list:
            wl, intensity, meta = load_spectrum(f)
            all_wl.append(wl)
            all_int.append(intensity)

        # 检查波长轴一致性
        ref_wl = all_wl[0]
        for wl in all_wl[1:]:
            if len(wl) != len(ref_wl) or np.max(np.abs(wl - ref_wl)) > 0.01:
                print("⚠ 波长轴不一致，改用逐文件拟合模式")
                return batch_fit(file_list, average_first=False, **kwargs)

        avg_intensity = np.mean(all_int, axis=0)
        std_intensity = np.std(all_int, axis=0)

        print(f"\n{'='*60}")
        print(f"  多次平均模式: {len(file_list)} 个光谱已平均")
        print(f"  平均后噪声降低约 {np.sqrt(len(file_list)):.1f} 倍")
        print(f"{'='*60}")

        result = fit_ruby(ref_wl, avg_intensity,
                         title=f'Averaged ({len(file_list)} spectra)', **kwargs)
        result['n_spectra'] = len(file_list)
        result['noise_reduction'] = np.sqrt(len(file_list))
        return result

    else:
        # 逐文件拟合
        all_results = []
        for f in file_list:
            print(f"\n--- {f.name} ---")
            try:
                wl, intensity, meta = load_spectrum(f)
                r = fit_ruby(wl, intensity, title=f.name, **kwargs)
                r['filename'] = f.name
                all_results.append(r)
            except Exception as e:
                print(f"  ✗ 拟合失败: {e}")

        if not all_results:
            raise RuntimeError("所有文件拟合失败")

        # 汇总统计
        r1_vals = [r['R1_nm'] for r in all_results]
        r2_vals = [r['R2_nm'] for r in all_results]
        fwhm_vals = [r['FWHM_R1'] for r in all_results]
        p_vals = [r['P_GPa'] for r in all_results]

        summary = {
            'n_spectra': len(all_results),
            'R1_mean': np.mean(r1_vals), 'R1_std': np.std(r1_vals),
            'R2_mean': np.mean(r2_vals), 'R2_std': np.std(r2_vals),
            'FWHM_R1_mean': np.mean(fwhm_vals), 'FWHM_R1_std': np.std(fwhm_vals),
            'P_mean': np.mean(p_vals), 'P_std': np.std(p_vals),
            'P_sem': np.std(p_vals) / np.sqrt(len(p_vals)),
            'individual': all_results,
            'models_used': list(set(r['model_label'] for r in all_results)),
        }

        print(f"\n{'='*60}")
        print(f"  批量统计 ({len(all_results)} 个文件)")
        print(f"  R1 峰位: {summary['R1_mean']:.4f} ± {summary['R1_std']:.4f} nm")
        print(f"  FWHM R1: {summary['FWHM_R1_mean']:.4f} ± {summary['FWHM_R1_std']:.4f} nm")
        print(f"  压  力 : {summary['P_mean']:.3f} ± {summary['P_std']:.3f} GPa "
              f"(SEM: {summary['P_sem']:.3f})")
        print(f"  采用模型: {', '.join(summary['models_used'])}")
        print(f"{'='*60}")

        return summary


# ╔══════════════════════════════════════════════════════════════════════╗
# ║  第八部分：使用示例                                                  ║
# ╚══════════════════════════════════════════════════════════════════════╝

if __name__ == '__main__':
    import glob

    DATA_DIR = Path(__file__).parent

    print("\n" + "▓" * 60)
    print("  红宝石压力标定 —— 最优拟合算法演示")
    print("▓" * 60)

    # ── 演示1: 单文件拟合（自动模型选择）──
    test_files = sorted(DATA_DIR.glob("Sample*1.9mW*1800*.csv"))
    if test_files:
        print(f"\n\n{'━'*60}")
        print("  演示1: 单文件拟合（自动模型选择）")
        print(f"{'━'*60}")
        wl, intensity, meta = load_spectrum(test_files[0])
        result = fit_ruby(wl, intensity, model='auto', baseline_order=1,
                         weighted=True, bootstrap_n=300,
                         title=test_files[0].name)

    # ── 演示2: 多次测量平均后拟合 ──
    if len(test_files) > 1:
        print(f"\n\n{'━'*60}")
        print("  演示2: 多次测量平均后拟合")
        print(f"{'━'*60}")
        avg_result = batch_fit(test_files, average_first=True,
                               model='auto', bootstrap_n=300)

    # ── 演示3: 逐文件拟合统计 ──
    test_files_1200 = sorted(DATA_DIR.glob("Sample*1.9mW*1200*.csv"))
    if test_files_1200:
        print(f"\n\n{'━'*60}")
        print("  演示3: 逐文件拟合统计 (1.9mW-1200)")
        print(f"{'━'*60}")
        batch_result = batch_fit(test_files_1200, average_first=False,
                                 model='auto', bootstrap_n=200)

    # ── 演示4: 高功率样本 ──
    high_power = sorted(DATA_DIR.glob("Sample*34mW*.csv"))
    if high_power:
        print(f"\n\n{'━'*60}")
        print("  演示4: 高功率样本")
        print(f"{'━'*60}")
        for f in high_power:
            wl, intensity, meta = load_spectrum(f)
            fit_ruby(wl, intensity, model='auto', baseline_order=2,
                    weighted=True, bootstrap_n=200, title=f.name)
