# -*- coding: utf-8 -*-
# @Time    : 2025/1/4 21:51
# @Author  : ljc
# @FileName: ulyss_fit_IDL_method.py
# @Software: PyCharm


# 1. 简介
"""
 Python conversion of the IDL uly_fit.pro .
目的:
    用 scipy 库中的数值优化方法, 推断 LASP 的恒星参数以及参数误差.
函数:
    1) uly_makeparinfo
    2) uly_fit_pparse
    3) uly_fitfunc_init
    4) uly_fitfunc
    5) uly_fit
解释:
    1) uly_makeparinfo 函数: 返回 cmp 字典结构中设置的待测光谱初始值等信息, 通过 uly_fit 函数调用.
    2) uly_fit_pparse 函数: 返回更新后的 losvd 参数 (降低 TGM 光谱分辨率)、以及在 cmp 结构中更新待
       测参数与误差, 通过 uly_fitfunc 与 uly_fit 函数调用.
    3) uly_fitfunc_init 函数: 返回初始化勒让德多项式参数 (降低 TGM 光谱与待测光谱的形状差异), 通过
       uly_fit 函数调用.
    4) uly_fitfunc 函数: 返回 TGM 光谱与待测光谱的流量差值 (即待优化的卡方函数, 使用 scipy 中的数
       值优化计算其极小值点), 通过 uly_fit 函数调用.
       注意: 
       4.1) LASP 设置的流量误差为 1, 因此并不是严格意义的卡方值.
       4.2) Python 版与 IDL 版一致, 对流量误差进行无偏估计, 并使用误差传播公式改正参数误差值.
       4.3) LASP 使用重复观测的统计方法估计参数误差.
    5) uly_fit 函数: 使用数值优化方法迭代更新 uly_fitfunc 获取最佳的恒星参数.
       5.1) 首先依次调用 uly_makeparinfo、uly_fitfunc_init 函数返回 losvd 参数初始值、恒星参数初始
       值、勒让德多项式系数初始值.
       5.2) 然后使用 scipy 中的数值优化迭代更新 uly_fitfunc 获取最佳的预测参数, 每次迭代时调用 uly_fit_pparse
       更新 losvd 参数与 cmp.
       注意: 
       5.3) 如果设置迭代剔除异常流量残差点, 则设置 uly_fit 函数中的 clean 参数为 True.
       5.4) 默认情况下 clean=False, 不迭代剔除异常流量残差点.
"""


# 2. 调库
import numpy as np
from scipy.special import eval_legendre
from uly_fit.uly_fit_lin import uly_fit_lin
from uly_read_lms.uly_spect_get import uly_spect_get
import time
from scipy.optimize import curve_fit
from uly_fit.robust_sigma import robust_sigma
import warnings
warnings.filterwarnings("ignore")


# 3. 获取 cmp 结构中设置的待测光谱初始值等信息
def uly_makeparinfo(cmp) -> list:

    """
        获取 cmp 结构中设置的待测光谱初始值等信息.

        输入参数:
        -----------
        cmp:
            为光谱字典结构, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.

        输出参数:
        -----------
        pinf:
             cmp 中设置的待测恒星大气物理参数信息.
    """

    # 3.1 获取 para 中的元素数量, 包含 Teff、log g、[Fe/H] 的设置, 因此长度为 3
    n_par = len(cmp['para']) if 'para' in cmp else 0
    if n_par == 0:
        raise ValueError("para is not specified!")

    # 3.2 创建恒星参数信息列表
    pinf = [{'value': 0.0, 
             'step': 1e-2,
             'limits': [0.0, 0.0], 
             'limited': [1, 1],
             'fixed': [0, 0]} for _ in range(n_par)
             ]

    # 3.3 填充恒星参数信息
    for j in range(n_par):
        pinf[j]['value'] = cmp['para'][j]['guess']         # 获取猜测值
        pinf[j]['step'] = cmp['para'][j]['step']           # 获取步长
        pinf[j]['limits'] = cmp['para'][j]['limits']       # 获取限制范围
        pinf[j]['limited'] = cmp['para'][j]['limited']     # 获取限制标志
        pinf[j]['fixed'] = cmp['para'][j]['fixed']         # 获取固定标志

    # 3.4 返回参数信息列表
    return pinf


# 4. 将 pars 的内容解析到 par_losvd 和 cmp 中
def uly_fit_pparse(pars=None, par_losvd=None, cmp=None, kmoment=None, error=None) -> tuple[list, dict]:

    """
        将 pars 的内容解析到 par_losvd 和 cmp 中.

        输入参数:
        -----------
        pars:
             参数数组, 前 2 个为 losvd 参数、第 3-5 个为待测恒星参数.
             LASP 中的 losvd 参数包括: cz, sigma (单位为像素), 用于降低 TGM 光谱分辨率以及计算 Rv.
        par_losvd:
                  losvd 参数.
        cmp:
            TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
        kmoment:
                losvd 的矩参数个数, LASP 设置为 2,  即 cz 与 sigma.
        error:
              待测参数的误差数组.

        输出参数:
        -----------
        par_losvd:
                  更新后的 losvd 参数, 用于降低 TGM 光谱分辨率.
        cmp:
            更新待测恒星参数、参数误差后的 cmp 字典.
    """

    # 4.1 如果 kmoment > 0 且 pars 不为空, 则将 pars 的前 kmoment 个元素赋值给 par_losvd
    if (kmoment > 0) and (pars is not None):
        par_losvd = pars[0: kmoment]

    # 4.2 将恒星参数存储到 cmp 结构中, 从 pars 的第 3 个元素开始 (因为前 2 个是 losvd 参数, 第 3-5 个是恒星参数)
    # 注意: 
    # 1) 我们尽可能保持与 IDL 一致的代码, 详情可参考 IDL 代码, 代价是代码可能冗余
    # 2) 后续版本将优化代码, 减少冗余, 如下述代码块可写成 cmp['para'][i]['value'] = pars[kmoment: kmoment + len(cmp['para'])][i]
    i0 = kmoment
    il = i0 + len(cmp['para']) - 1    # 2+3-1=4
    if il >= i0:
        for i in range(len(cmp["para"])):
            cmp['para'][i]['value'] = pars[i0: il + 1][i]

    # 4.3 将恒星参数误差存储到 cmp 结构中, 这里仅在迭代剔除异常流量残差时使用
    if error is not None:
        if il >= i0:
            for i in range(len(cmp["para"])):
                cmp['para'][i]['error'] = error[i0: il + 1][i]

    # 4.4 返回 losvd 参数、以及更新恒星参数与参数误差后的 cmp
    return par_losvd, cmp


# 5. 初始化用于 uly_fit_lin 的缓存信息
def uly_fitfunc_init(spec=None, mpoly=None, cmp=None, mdegree=None, modecvg=None) -> tuple[dict, dict, int]:

    """
        初始化用于 uly_fit_lin 的缓存信息.

        输入参数:
        -----------
        spec:
            待测光谱字典结构.
        mpoly:
              字典结构. 包含: {lmdegree, mpolcoefs, poly, leg_array}.
              1) lmdegree: 勒让德多项式的最大阶数 (即 n, LAMOST 使用 50 阶).
                 1.1) 迭代过程中保持不变 (伪连续谱系数为正. 或负数时推断值为 -9999).
                 1.2) 迭代过程中可减少 (伪谱系数为负数, 减少勒让德多项式阶数使得伪谱系数为正).
              2) mpolcoefs: 勒让德多项式的系数, 即 (c0, c1, c2, ..., cn), 形状是 (lmdegree+1, 1), 该系数使用最小二乘法求解.
              3) leg_array: 勒让德多项式关于 x (可视为 lambda_i) 的值组成的矩阵, 形状是 (流量维度, lmdegree+1), 各元素为 (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)).
                 注意: 3.1) lambda_i 是第 i 个像素的波长;
                      3.2) 为了便于计算, LASP-MPFit 或 LASP-Adam-GPU 中使用等间隔 [-1, 1) 中的序列表示光谱波长, 而不是真实波长点;
                      3.3) Pn(lambda_i) 是第 n 阶勒让德多项式在 lambda_i 处的值.
              4) poly: 各波长处的多项式值组成的伪谱 (流量改正因子), 即 B(lambda_i), 形状是 (流量维度, 1).
                 注意: 4.1) n 阶勒让德多项式建立的伪谱 (流量改正因子) 中的第 i 个像素值: B(lambda_i) = c0 + c1 * P1(lambda_i) + c2 * P2(lambda_i) + ... + cn * Pn(lambda_i).
        cmp:
            TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
        mdegree:
                勒让德多项式阶数.
                注意: 
                1) LASP 使用的 50 阶勒让德多项式.
                2) 默认情况下 mdegree=10.
        modecvg:
                指定拟合方法的收敛模式. LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码.
                注意:
                1) modecvg=0 (默认选项), 这是最快的, 但如果解错失, 问题可能发生.
                2) modecvg=1 (每次迭代仅完成一次), 为每个 LM 迭代计算导数.
                3) modecvg=2, uly_fit_lin 始终收敛, 但速度较慢.

        输出参数:
        -----------
        mpoly:
              初始化勒让德多项式相关系数.
        cmp:
            传递 cmp 字典结构信息.
        modecvg:
                收敛模式.
    """

    # 5.1 检查输入数据是否完整
    if cmp is None:
        raise ValueError('CMP dictionary is not specified!')
    if spec is None:
        raise ValueError('SPECTRUM parameter must be specified!')
    if mdegree is None:
        raise ValueError('MDEGREE parameter must be specified!')

    # 5.2 获取待测光谱像素数点
    npix = spec['data'].shape[0]

    # 5.3 初始化/重新初始化 mpoly
    # 5.3.1 如果 mpoly 字典结构体不为空, 则检查 mpoly 字典结构体中的勒让德多项式最大次数与勒让德多项式值的形状是否与 mdegree 和 npix 一致
    if mpoly is not None:
        if (mpoly['lmdegree'] != mdegree) | (mpoly['poly'].shape[0] != npix):
            mpoly = {}
    # 5.3.2 如果 mpoly 字典结构体为空, 则创建 mpoly 字典
    if mpoly is None:
        mpoly = {'lmdegree': mdegree,                           # 勒让德多项式的最大次数
                 'mpolcoefs': np.zeros(mdegree+1),              # 勒让德多项式的系数
                 'poly': np.ones(shape=npix),                   # 伪连续谱
                 'leg_array': np.ones(shape=(npix, mdegree+1))  # 勒让德多项式关于 x 的值
                 }
        # 5.3.2.1 计算勒让德多项式
        # 1) 设置勒让德多项式自变量范围为 [-1, 1]
        # 2) 代入勒让德多项式, 计算勒让德多项式关于自变量 x 的值
        x = 2.0 * np.arange(npix, dtype=np.float64) / npix - 1.0
        mpoly['leg_array'] = eval_legendre(np.arange(mdegree+1), x[:, np.newaxis])

    # 5.3.2.2 初始化 mpoly 的其他属性
    mpoly["mpolcoefs"][0] = 1    # 勒让德多项式的常数项系数为 1

    # 5.3.3 设置收敛模式, LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码
    if modecvg is None:
        modecvg = 0

    # 5.3.4 返回初始的勒让德多项式值、传递 cmp 信息、以及一个不重要的参数 modecvg
    return mpoly, cmp, modecvg


# 6. 由 scipy 中的优化器迭代优化的流量误差函数
def uly_fitfunc(pars, kpen=None, polpen=None, adegree=None, kmoment=None, voff=None, 
                signalLog=None, goodpixels=None, lum_weight=None, outpixels=None, sampling_function=None,
                mpoly=None, cmp=None, modecvg=None, addcont=None, allow_polynomial_reduction=False) -> tuple[np.ndarray, np.ndarray, dict, dict]:

    """
       由 scipy 中的优化器迭代优化的流量误差函数.

       输入参数:
       -----------
       pars:
            待测的参数数组, 前 2 个是运动学高斯核参数 (降低分辨率), 第 3-5 是 3 个大气物理参数.
       kpen:
            此参数会将 (h3, h4, ...) 的测量值偏向零, 除非其包含项显著减少了流量拟合误差.
            注意: 
            1) 默认情况下, kpen=0, 表示未启用惩罚项 (LASP 中默认设置为 0).
            2) 如果设置为严格正值, 解 (包括 cz 和 sigma) 会减少噪声. 使用惩罚时, 建议使用蒙特卡洛模拟测试 kpen 的选择.
            3) kpen 的值范围应在 0.5 到 1.0 之间, 作为好的初始猜测.
       polpen:
              乘法多项式的偏置水平. 此关键词可用于减少乘法多项式中不重要项的影响.
              注意: 
              1) 默认情况下不应用偏置. 如果某些系数的绝对值小于 polpen 倍的统计
                 误差, 这些系数会通过因子 (abs(coef)/(polpen*err))^2 被抑制.
              2) 该功能仅在 mdegree>0 时有效, polpen=2 是一个合理的选择.
       adegree:
               用于修正 TGM 模板光谱形状的加法勒让德多项式的阶数. 
               注意:
               1) 在拟合过程中, 默认不使用任何加法多项式.
               2) 如果要禁用加法多项式, 请设置 adegree=-1, LASP 设置为 -1.
       kmoment:
               高斯-厄米特矩的阶数.
               设置为 2 时仅拟合 [cz, sigma].
               设置为 4 时同时拟合 [cz, sigma, h3, h4].
               设置为 6 时同时拟合 [cz, sigma, h3, h4, h5, h6].
               注意: 
               1) LASP 设置为 2.
       voff:
            模型光谱与待测光谱的波长偏移, 单位为 km/s.
       signalLog:
                 包含待测光谱的结构体 (见 CMP), 模型光谱和待测光谱的波长必须是 ln 对数化的, 相关标签在 uly_spect_alloc 文档中有描述.
       goodpixels:
                  好像素点的索引值.
       quiet:
             是否抑制屏幕上打印的消息.
             注意:
             1) quiet=True, 表示抑制屏幕上打印的消息.
             2) quiet=False, 表示不抑制屏幕上打印的消息.
       lum_weight:
                  计算流量权重及其误差. 如果设置, 计算 e_weight 和 l_weight.
       outpixels:
                 goodpixels.
       sampling_function:
                        插值方法. 可输入 "splinf", "cubic", "slinear", "quadratic", "linear". 默认使用 "linear" 插值方法.
       mpoly:
              字典结构. 包含: {lmdegree, mpolcoefs, poly, leg_array}.
              1) lmdegree: 勒让德多项式的最大阶数 (即 n, LAMOST 使用 50 阶).
                 1.1) 迭代过程中保持不变 (伪连续谱系数为正. 或负数时推断值为 -9999).
                 1.2) 迭代过程中可减少 (伪谱系数为负数, 减少勒让德多项式阶数使得伪谱系数为正).
              2) mpolcoefs: 勒让德多项式的系数, 即 (c0, c1, c2, ..., cn), 形状是 (lmdegree+1, 1), 该系数使用最小二乘法求解.
              3) leg_array: 勒让德多项式关于 x (可视为 lambda_i) 的值组成的矩阵, 形状是 (流量维度, lmdegree+1), 各元素为 (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)).
                 注意: 3.1) lambda_i 是第 i 个像素的波长;
                      3.2) 为了便于计算, LASP-MPFit 或 LASP-Adam-GPU 中使用等间隔 [-1, 1) 中的序列表示光谱波长, 而不是真实波长点;
                      3.3) Pn(lambda_i) 是第 n 阶勒让德多项式在 lambda_i 处的值.
              4) poly: 各波长处的多项式值组成的伪谱 (流量改正因子), 即 B(lambda_i), 形状是 (流量维度, 1).
                 注意: 4.1) n 阶勒让德多项式建立的伪谱 (流量改正因子) 中的第 i 个像素值: B(lambda_i) = c0 + c1 * P1(lambda_i) + c2 * P2(lambda_i) + ... + cn * Pn(lambda_i).
       cmp:
           TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
           cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
       modecvg:
               指定拟合方法的收敛模式. LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码.
               注意:
               1) modecvg=0 (默认选项), 这是最快的, 但如果解错失, 问题可能发生.
               2) modecvg=1 (每次迭代仅完成一次), 为每个 LM 迭代计算导数.
               3) modecvg=2, uly_fit_lin 始终收敛, 但速度较慢.
       addcont:
               加法连续体的数组.
               注意:
               1) 在拟合过程中, 默认不使用任何加法多项式.
       allow_polynomial_reduction:
                                  是否允许多项式阶数减少. 默认值为 False, 即不允许多项式阶数减少.
                                  如果设置为 True, 则允许多项式阶数减少.
                                  注意: 伪连续谱不应该小于 0, 如果小于 0, 提供两种处理方法:
                                  1) 为避免伪连续谱存在负值, 循环减少多项式阶数, LASP IDL 版本采用这种方法.
                                  2) 直接认为该光谱质量较差, 参数推断失败, Python 版本默认采用这种方法 (因为连续谱为负值, 大概率由于光谱流量存在负值).

        输出参数：
        -----------
        err:
            TGM 模型光谱与待测光谱的流量差值.
        bestfit:
                TGM 光谱.
        mpoly:
               乘法勒让德多项式结构体.
        cmp:
            由 uly_fit_pparse 函数更新后的 cmp 结构.
    """

    # 6.1 检查输入参数是否包含 NaN
    if pars is not None:
        if np.any(np.isnan(pars)):
            raise ValueError('(pars) array contains NaN!')

    # 6.2 如果未提供 outpixels, 则使用 goodpixels
    if outpixels is None:
        outpixels = goodpixels

    # 6.3 根据输入参数 pars, 获取 losvd 参数、待测恒星参数
    # 1) pars 的前 2 个为 losvd 参数
    # 2) pars 的第 3-5 个为恒星参数
    par_losvd, cmp = uly_fit_pparse(pars=pars,         # 待测参数
                                    par_losvd=None,    # losvd 参数
                                    cmp=cmp,           # cmp 组件字典
                                    kmoment=kmoment    # losvd 参数数量
                                    )

    # 6.4 调用 uly_fit_lin 函数, 获取指定参数下的 TGM 光谱
    # 1) 首先在 uly_fit_lin 函数中调用 uly_tgm_eval 函数获取与待测光谱相同流量点的 TGM 光谱 (使用 uly_spect_logrebin 函数插值)
    # 2) 然后在 uly_fit_lin 函数中调用 convol 函数获取指定 par_losvd 参数下的低分辨率 TGM 光谱
    # 3) 然后在 uly_fit_lin 函数中调用 uly_fit_lin_weight 函数获取加权后的 TGM 光谱
    # 4) 最后在 uly_fit_lin 函数中调用 uly_fit_lin_mulpol 函数获取乘以勒让德多项式 (mpoly) 的 TGM 光谱
    bestfit, mpoly = uly_fit_lin(adegree=adegree,                                              # 加法勒让德多项式阶数
                                 polpen=polpen,                                                # 乘法勒让德多项式偏置水平
                                 voff=voff,                                                    # 速度偏移
                                 par_losvd=par_losvd,                                          # losvd 参数
                                 cmp=cmp,                                                      # cmp 字典
                                 goodPixels=goodpixels,                                        # 有效像素列表
                                 sampling_function=sampling_function,                          # 插值方法
                                 SignalLog=signalLog,                                          # 待测光谱字典结构
                                 mpoly=mpoly,                                                  # 勒让德多项式字典结构
                                 addcont=addcont,                                              # 加法多项式数组
                                 modecvg=modecvg,                                              # 模式收敛, LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码
                                 lum_weight=lum_weight,                                        # 计算流量权重及其误差
                                 allow_polynomial_reduction=allow_polynomial_reduction         # 是否允许多项式阶数减少
                                 )
    bestfit = bestfit.reshape(-1)

    # 6.5 检查模式收敛, LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码
    if modecvg == 1:
        modecvg = 3

    # 6.6 计算模型和观测之间的加权偏差 (卡方值)
    # 注意:
    # 1) LASP 没有使用待测光谱的流量误差, 因此 signalLog['err'] 为全 1 数组
    # 2) 因此基于卡方优化得到的参数误差需要改正
    err = (signalLog['data'][outpixels] - bestfit[outpixels]) / signalLog['err'][outpixels]

    # 6.7 如果 kmoment>2 且 kpen!=0, 则对流量进行惩罚
    if kmoment > 2 and kpen != 0:
        sigma = robust_sigma(err, zero=True)
        err = err + kpen * sigma * np.sqrt(np.sum(pars[2:kmoment] ** 2))

    # 6.8 返回 TGM 光谱与待测光谱的流量差值、TGM 光谱、勒让德多项式字典结构、cmp 字典结构
    return err, bestfit, mpoly, cmp


# 7. 由 scipy 中的数值优化方法推断待测的 2 个 losvd 参数、3 个恒星大气物理参数、以及它们的误差
def uly_fit_IDL_method(signalLog, cmp=None, kmoment=None, kguess=None, kfix=None, klim=None,
                       kpen=None, adegree=None, mdegree=None, sampling_function=None, polpen=None, clean=False,
                       modecvg=None, allow_polynomial_reduction=False, quiet=False, 
                       full_output=False, plot_fitting=False) -> tuple[float, float, float, float, float, float, float, float, float, float]:

    """
       在 IDL 中, 由 ulyss 直接调用的函数. 该函数组合 uly_makeparinfo、uly_fit_pparse、uly_fitfunc_init、uly_fitfunc, 
       使用 scipy 中的数值优化方法推断待测的 2 个 losvd 参数、3 个恒星大气物理参数、以及它们的误差.

       输入参数:
       -----------
       signalLog:
                 包含待测光谱的结构体 (见 CMP), 模型光谱和待测光谱的波长必须是 ln 对数化的, 相关标签在 uly_spect_alloc 文档中有描述.
       cmp:
           TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
           cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
       kmoment:
               高斯-厄米特矩的阶数.
               设置为 2 时仅拟合 [cz, sigma].
               设置为 4 时同时拟合 [cz, sigma, h3, h4].
               设置为 6 时同时拟合 [cz, sigma, h3, h4, h5, h6].
               注意: 
               1) LASP 设置为 2.
       kguess:
              待测参数的初始值.
       kfix:
            是否固定待测 losvd 参数. 1 表示在最小化过程中固定对应的 losvd 参数. 例如, 若要固定速度离散度, 可以指定 kfix=[0,1].
       klim:
            待测参数的上下边界.
       kpen:
            此参数会将 (h3, h4, ...) 的测量值偏向零, 除非其包含项显著减少了流量拟合误差.
            注意: 
            1) 默认情况下, kpen=0, 表示未启用惩罚项 (LASP 中默认设置为 0).
            2) 如果设置为严格正值, 解 (包括 cz 和 sigma) 会减少噪声. 使用惩罚时, 建议使用蒙特卡洛模拟测试 kpen 的选择.
            3) kpen 的值范围应在 0.5 到 1.0 之间, 作为好的初始猜测.
       adegree:
               用于修正 TGM 模板光谱形状的加法勒让德多项式的阶数. 
               注意:
               1) 在拟合过程中, 默认不使用任何加法多项式.
               2) 如果要禁用加法多项式, 请设置 adegree=-1, LASP 设置为 -1.
       mdegree:
               勒让德多项式阶数.
               注意: 
               1) ULySS 默认是 10.
               2) LASP 设置为 50.
               3) 不同任务, mdegree 的设置不同, 可参考 best_mdegree.ipynb 文件对 mdegree 进行设置.
       polpen:
              乘法多项式的偏置水平. 此关键词可用于减少乘法多项式中不重要项的影响.
              注意: 
              1) 默认情况下不应用偏置. 如果某些系数的绝对值小于 polpen 倍的统计
                误差, 这些系数会通过因子 (abs(coef)/(polpen*err))^2 被抑制.
              2) 该功能仅在 mdegree>0 时有效, polpen=2 是一个合理的选择.
       clean:
             是否以迭代方式检测并剪切 TGM 光谱与待测光谱流量残差的离群值.
             注意:
             1) Clean=True, 表示以迭代方式检测并剪切 TGM 光谱与待测光谱流量残差的离群值.
             2) Clean=False, 表示不进行离群值检测与剪切.
       modecvg:
               指定拟合方法的收敛模式. LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码.
               注意:
               1) modecvg=0 (默认选项), 这是最快的, 但如果解错失, 问题可能发生.
               2) modecvg=1 (每次迭代仅完成一次), 为每个 LM 迭代计算导数.
               3) modecvg=2, uly_fit_lin 始终收敛, 但速度较慢.
       allow_polynomial_reduction:
                                  是否允许多项式阶数减少. 默认值为 False, 即不允许多项式阶数减少.
                                  如果设置为 True, 则允许多项式阶数减少.
                                  注意: 伪连续谱不应该小于 0, 如果小于 0, 提供两种处理方法:
                                  1) 为避免伪连续谱存在负值, 循环减少多项式阶数, LASP IDL 版本采用这种方法.
                                  2) 直接认为该光谱质量较差, 参数推断失败, Python 版本默认采用这种方法 (因为连续谱为负值, 大概率由于光谱流量存在负值).
       quiet:
             是否抑制屏幕上打印的消息.
             注意:
             1) quiet=True, 表示抑制屏幕上打印的消息.
             2) quiet=False, 表示不抑制屏幕上打印的消息.
       full_output:
                   是否返回所有参数的拟合结果, 以及拟合信息.
                   注意:
                   1) full_output=True, 表示返回所有参数的拟合结果, 以及拟合信息.
                   2) full_output=False, 表示仅返回 Rv、Teff、log g、[Fe/H] 以及误差推断值.
       plot_fitting:
                   是否绘制光谱的拟合流量残差图.
                   注意:
                   1) plot_fitting=True, 表示绘制拟合流量残差图.
                   2) plot_fitting=False, 表示不绘制拟合流量残差图.

        输出参数：
        -----------
        Rv, Teff, logg, FeH, Rv_err, Teff_err, logg_err, FeH_err, used_time, loss
        注意:
        1) used_time 表示推断 1 条光谱的参数所用的时间, 单位为秒.
        2) loss 表示拟合的流量残差的均方根误差.
    """

    # 7.1 初始信息
    # 7.1.1 定义光速常数 (km/s), 从 signalLog 字典中获取待测光谱, 获取待测光谱的像素点数量, 检查是否设置了光谱流量误差
    c0, galaxy, npix, have_noise  = 299792.458, signalLog["data"], signalLog["data"].shape[0], 1
    # 7.1.2 检查是否设置了光谱流量误差, 如果没有设置, 则创建全 1 数组作为光谱流量误差
    if (len(signalLog["err"]) == 0) | (signalLog["err"] == ""):
        have_noise = 0
        signalLog["err"] = np.ones_like(galaxy)
    # 7.1.3 获取光谱流量误差
    noise = signalLog["err"]

    # 7.2 获取 msk 列表
    # 1) 1 是好像素点, 0 是坏像素点
    # 2) uly_spect_get 返回 WAVERANGE, GOODPIX, HDR, MASK, 因此使用第 3 个位置获取 msk 信息
    msk = uly_spect_get(signalLog,  # 输入光谱字典结构
                        MASK=True   # 获取 msk 信息
                        )[3]
    # 7.2.1 如果设置了 cmp 字典, 则将 cmp 字典中的 mask 信息与 msk 相乘
    if cmp is not None:
        if ("mask" in cmp) and ((cmp["mask"] is not None) or (cmp["mask"] != "")):
            if len(cmp["mask"]) > 0:
                msk *= cmp["mask"]
    # 7.2.2 获取好像素的索引值
    goodpix = uly_spect_get(signalLog,     # 输入光谱字典结构
                            GOODPIX=["1"]  # 获取好像素列表
                            )[1]

    # 7.3 检查并处理 NaN 值
    # 7.3.1 检查光谱中的 NaN 值
    nan_indices = np.where(~np.isfinite(galaxy[goodpix]))[0].tolist()
    if len(nan_indices) > 0:
        if quiet is True:
            print(f'{len(nan_indices)} 个 NaN 值在光谱中被屏蔽!')
        msk[goodpix[nan_indices], :] = 0
    # 7.3.2 检查误差光谱中的 NaN 值
    nan_indices = np.where(~np.isfinite(noise[goodpix]))[0].tolist()
    if len(nan_indices) > 0:
        if quiet is True:
            print(f'{len(nan_indices)} 个 NaN 值在误差光谱中被屏蔽!')
        msk[goodpix[nan_indices], :] = 0
    # 7.3.3 获取所有好的像素点
    goodPixels0 = np.where(msk == 1)[0].tolist()
    if len(goodPixels0) == 0:
        raise ValueError('No good pixels left!')

    # 7.4 计算速度偏移
    # 7.4.1 计算速度偏移
    voff = (cmp["start"] - signalLog["start"]) * c0
    # 7.4.2 计算速度刻度
    velScale = c0 * signalLog["step"]

    # 7.5 加型、乘法多项式参数初始化
    # 7.5.1 如果 adegree 为 None, 则设置为 -1, 表示不使用加型多项式
    if adegree is None:
        adegree = -1
    # 7.5.2 如果 adegree 不为 None, 则使用 adegree, LASP 使用 -1
    else:
        adegree = adegree
    # 7.5.3 如果 mdegree 为 None, 则设置为 10, 表示使用 10 阶勒让德多项式
    if mdegree is None:
        mdegree = 10
    # 7.5.4 如果 mdegree 不为 None, 则使用 mdegree, LASP 使用 50
    else:
        mdegree = mdegree

    # 7.6 检查像素范围是否有效
    if np.max(goodPixels0) > len(galaxy) - 1:
        raise ValueError('goodpixels out of range!')

    # 7.7 设置 kpen 默认值
    if kpen is None:
        kpen = 0.7 * np.sqrt(500 / len(goodPixels0))

    # 7.8 检查 kmoment 参数
    # 7.8.1 如果 kmoment 为 None, 则设置为 2, 表示拟合 [cz, sigma]
    if kmoment is None:
        kmoment = 2
    # 7.8.2 如果 kmoment 不为 None, 则检查 kmoment 是否在 [0, 2, 4, 6] 中
    elif kmoment not in [0, 2, 4, 6]:
        raise ValueError('KMOMENT should be 0, 2, 4 or 6!')

    # 7.9 检查 kguess 参数
    # 7.9.1 如果 kmoment >= 2 且 kguess 为 None 或 kguess 的长度小于 2, 则抛出异常
    if (kmoment >= 2) and (kguess is None or len(kguess) < 2):
        raise ValueError('kguess must have two elements [V, sigma]!')

    # 7.10 设置 kfix 参数
    # 7.10.1 如果 kfix 为 None, 则设置为全 0 的数组
    nst = len(kguess)
    if kfix is None:
        kfix = np.zeros(nst, dtype=int)
    # 7.10.2 如果 kfix 的长度小于 kguess 的长度, 则将 kfix 的长度扩展为 kguess 的长度
    elif len(kfix) < nst:
        kfix = np.concatenate([kfix, np.zeros(nst - len(kfix), dtype=int)])
    # 7.10.3 如果 kfix 的长度大于 kguess 的长度, 则将 kfix 的长度截断为 kguess 的长度
    elif len(kfix) > nst:
        kfix = kfix[: nst]

    # 7.11 losvd 参数的限制设置
    # 7.11.1 如果 kmoment > 0, 则设置 klimits
    if kmoment > 0:
        klimits = np.zeros((2, kmoment))
        # 7.11.2 设置速度限制
        klimits[0:2, 0] = (kguess[0] + np.array([-2e3, 2e3]))  # 第 0 行, 前 2 列
        # 7.11.3 设置 sigma 限制
        if kmoment >= 2:
            klimits[0:2, 1] = [0.3 * velScale, 1e3]  # 第 1 行, 前 2 列
        # 7.11.4 设置 h3, h4 等的限制
        for k in range(2, kmoment):
            klimits[0:2, k] = [-0.3, 0.3]
        # 7.11.5 如果给定了实际限制, 则覆盖默认值
        if klim is not None:
            klimits[0] = klim
        # 7.11.6 诊断超出边界的情况
        for k in range(kmoment):
            if (kfix[k] == 0) and (kguess[k] < klimits[0, k]):
                raise ValueError(f'Guess on kinematic moment {k}: {kguess[k]} '
                                 f'is lower than the limit: {klimits[0, k]}')
            elif kfix[k] == 0 and kguess[k] > klimits[1, k]:
                raise ValueError(f'Guess on kinematic moment {k}: {kguess[k]} '
                                 f'is higher than the limit: {klimits[1, k]}')
        # 7.11.7 将 klimits 转换为像素单位
        klimits[0: 2, 0] = np.log(1 + klimits[0: 2, 0] / c0) / signalLog["step"]
        if kmoment >= 2:
            klimits[0: 2, 1] = np.log(1 + klimits[0: 2, 1] / c0) / signalLog["step"]

    # 7.12 准备恒星参数信息
    # 注意:
    # 1) Python 版支持 1 次初始化 1 组参数, 即使用 CFI 初始值作为迭代起始点, 而不是使用多组初始值
    # 2) 恒星参数信息包括 Rv、Teff、log g、[Fe/H]
    # 7.12.1 获取好像素点
    goodPixels = goodPixels0.copy()
    # 7.12.2 测试参数猜测值是否在允许的范围内
    for k in range(len(cmp["para"])):
        # 7.12.2.1 检查下限
        if cmp["para"][k]["guess"] < cmp["para"][k]["limits"][0]:
            raise ValueError(f'guess on {cmp["para"][k]["name"]} : '
                             f'{cmp["para"][k]["guess"]} '
                             f'is below the limit : {cmp["para"][k]["limits"][0]}')
        # 7.12.2.2 检查上限
        elif cmp["para"][k]["guess"] > cmp["para"][k]["limits"][1]:
            raise ValueError(f'guess on {cmp["para"][k]["name"]} : '
                             f'{cmp["para"][k]["guess"]} '
                             f'is above the limit : {cmp["para"][k]["limits"][1]}')
    # 7.12.3 创建 parinfok: 包含 losvd 的 parinfo 部分
    if kmoment == 0:
        parinfok = None
    if kmoment > 0:
        # 7.12.3.1 创建 parinfok 字典结构
        parinfok = [{'value': 0.0, 'step': 1e-2, 'limits': [0.0, 0.0], 'limited': [1, 1], 'fixed': 0} for _ in range(kmoment)]
        # 7.12.3.2 设置第一个参数(速度)
        # 7.12.3.2.1 转换速度为像素
        parinfok[0]['value'] = np.log(1 + kguess[0] / c0) / signalLog["step"]
        # 7.12.3.2.2 设置速度限制
        parinfok[0]['limits'] = klimits[:, 0]
        # 7.12.3.2.3 如果径向速度固定, 则设置固定
        if kfix[0] == 1:
            parinfok[0]['fixed'] = 1
            parinfok[0]['step'] = 0
            parinfok[0]['limits'] = [parinfok[0]['value'], parinfok[0]['value']]
    # 7.12.4 设置第二个运动学参数(如果存在)
    if kmoment > 1:
        # 7.12.4.1 转换速度为像素
        parinfok[1]['value'] = np.log(1 + kguess[1] / c0) / signalLog["step"]
        # 7.12.4.2 设置速度限制
        parinfok[1]['limits'] = klimits[:, 1]
        # 7.12.4.3 如果速度弥散固定, 则设置固定
        if kfix[1] == 1:
            parinfok[1]['fixed'] = 1
            parinfok[1]['step'] = 0
            parinfok[1]['limits'] = [parinfok[1]['value'], parinfok[1]['value']]
        # 7.12.4.4 检查波长范围
        if cmp["npix"] <= ((abs(voff) + abs(kguess[0]) + 5 * kguess[1]) / velScale):
            raise ValueError('Wavelength range is too small, or velocity shift too big!')
    # 7.12.5 设置更高阶的运动学参数(如果存在)
    if kmoment > 2:
        # 7.12.5.1 设置 h3, h4 等参数
        for k in range(2, kmoment):
            parinfok[k]['value'] = kguess[k]
            # 7.12.5.1.1 设置 h3, h4 等参数限制
            parinfok[k]['limits'] = [-0.3, 0.3]
            # 7.12.5.1.2 设置 h3, h4 等参数步长
            parinfok[k]['step'] = 1e-3
        # 7.12.5.2 处理固定参数
        fixh = np.where(kfix == 1)[0].tolist()
        if len(fixh) > 0:
            for idx in fixh:
                parinfok[idx]['fixed'] = 1
                parinfok[idx]['step'] = 0
                parinfok[idx]['limits'] = [kguess[idx], kguess[idx]]

    # 7.13 创建组件部分的 parinfo, 包含恒星参数信息
    # 7.13.1 获取恒星大气物理参数信息
    parinfo = uly_makeparinfo(cmp)
    # 7.13.2 合并 parinfo 和 parinfok
    if parinfok is not None:
        if parinfo is not None:
            # 7.13.2.1 合并 parinfo 和 parinfok, 更新 parinfo
            # 其中 parinfok 包含 losvd 的参数信息
            parinfo = parinfok + parinfo
        else:
            # 7.13.2.2 如果 parinfo 为 None, 则直接使用 parinfok
            parinfo = parinfok

    # 7.14 以一组猜测值为起点计算函数极小值点
    # 注意:
    # 1) LASP-MPFit 以 CFI 初始值作为迭代起始点, 如果 CFI 没被提供, 则使用多组初始值
    # 2) LASP-CurveFit, LASP-Adam-GPU 以 CFI 初始值作为迭代起始点, 而不是使用多组初始值, 如果 CFI 没被提供, 则使用一组固定值作为初始值
    # 3) IDL 与 Python 均包含 clean 模式 与 No Clean 模式, 但 LASP-Adam-GPU 采用一致的策略, 同时 MASK 多条光谱
    # 4) No Clean 模式: 找到最小值后直接跳出循环, 默认不清理流量残差异常值
    # 5) Clean 模式: 找到最小值后, 清理流量残差异常值并重复最小化, 最多清理 11 次, 每次判断使用 3、4、5、7 倍标准差清理
    # 7.14.1 设置清理迭代次数, 默认不清理
    nclean = 0
    # 7.14.2 如果需要 clean 模式, 则设置清理迭代次数
    if clean is True:
        nclean = 10
    # 7.14.3 IDL 设置的优化容差在 scipy 的数值优化中不可取, 需要设置的再小一些
    # 注意:
    # 1) 默认不清理流量残差异常值, 即 No Clean 模式, 优化容差为 1e-5
    # 2) 如果需要清理流量残差异常值, 即 Clean 模式, 优化容差为 1e-2
    ftol = 1e-5 if nclean == 0 else 1e-2
    xtol = 1e-10 if nclean == 0 else 1e-8

    # 7.14.4 对每个清理迭代进行循环. 
    # 注意: 
    # 1) 我们尽可能保证与 IDL 版的一致性, 因此代码可能比较冗余, 后续更新应考虑代码优化
    # 7.14.4.1 记录推断出 1 条光谱最佳恒星参数的起始时间
    time0 = time.time()
    for j in range(nclean + 1):
        # 7.14.4.2 初始化
        mpoly, cmp, modecvg = uly_fitfunc_init(spec=signalLog,  # 输入光谱字典结构
                                               mpoly=None,      # 勒让德多项式字典结构
                                               cmp=cmp,         # 字典组件
                                               mdegree=mdegree, # 勒让德多项式阶数
                                               modecvg=modecvg  # 迭代模式, LASP 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码
                                               )
        # 7.14.4.3 初始化 [2 个 losvd 参数、3 个恒星大气物理参数]
        value = [parinfo[i]["value"] for i in range(len(parinfo))]

        # 7.14.4.4 构建 scipy 中的优化器待优化的目标函数
        def min_fun(temp, *pars) -> np.ndarray:
            # 7.14.4.4.1 调用 uly_fitfunc 计算流量残差, 由于 uly_fitfunc 返回值为 (resc0, bestfit, mpoly, cmp), 所以取第一个元素
            resc0 = uly_fitfunc(pars=pars,                                               # 待测参数
                                kpen=kpen,                                               # 惩罚因子
                                cmp=cmp,                                                 # 字典组件
                                adegree=adegree,                                         # 勒让德多项式阶数
                                signalLog=signalLog,                                     # 光谱字典数据
                                mpoly=mpoly,                                             # 勒让德多项式字典结构
                                goodpixels=goodPixels,                                   # 有效像素列表
                                sampling_function=sampling_function,                     # 插值方法
                                voff=voff / velScale,                                    # 速度偏移
                                kmoment=kmoment,                                         # 运动学参数数量
                                allow_polynomial_reduction=allow_polynomial_reduction    # 是否允许多项式阶数减少
                                )[0]
    
            # 7.14.4.4.2 返回流量残差
            return resc0

        # 7.14.4.5 计算待测参数 (自由非线性参数) 的数量
        nlin_free = 0
        for k in range(len(parinfo)):
            if parinfo[k]['fixed'] == 0:
                nlin_free += 1

        # 7.14.4.6 如果有待测参数, 则开始推断
        if nlin_free > 0:
            # 7.14.4.6.1 使用 scipy 中的优化器
            # 7.14.4.6.1.1 使用 least_squares 方法
            """
            # 使用 least_squares 方法
            result = least_squares(min_fun,
                                    x0=value,               # value 前两个是降低分辨率的, 第3-5个是大气参数, 后面几个勒让德多项式系数
                                    method='lm',            # Levenberg-Marquardt 算法
                                    ftol=1e-8, xtol=1e-8)   # 收敛标准
            res = result.x
            error = result.cov_x if hasattr(result, 'cov_x') else None
            bestnorm = result.fun
            mpfstat = result.status
            ncalls = result.nfev
            """

            # 7.14.4.6.2 使用 curve_fit, 返回参数推断均值与标准差
            # 7.14.4.6.2.1 设置待测参数边界值
            bounds = np.zeros(shape=(2, 5))
            for lim in range(len(parinfo)):
                bounds[0, :] = np.array([parinfo[0]["limits"][0],   # 设置 losvd 第 1 个参数下限 (核函数均值)
                                         parinfo[1]["limits"][0],   # 设置 losvd 第 2 个参数下限 (核函数标准差)
                                         parinfo[2]["limits"][0],   # 设置恒星大气物理参数第 1 个参数下限 (Teff)
                                         parinfo[3]["limits"][0],   # 设置恒星大气物理参数第 2 个参数下限 (logg)
                                         parinfo[4]["limits"][0]    # 设置恒星大气物理参数第 3 个参数下限 ([Fe/H])
                                         ])
                bounds[1, :] = np.array([parinfo[0]["limits"][1],   # 设置 losvd 第 1 个参数上限 (核函数均值)
                                         parinfo[1]["limits"][1],   # 设置 losvd 第 2 个参数上限 (核函数标准差)
                                         parinfo[2]["limits"][1],   # 设置恒星大气物理参数第 1 个参数上限 (Teff)
                                         parinfo[3]["limits"][1],   # 设置恒星大气物理参数第 2 个参数上限 (logg)
                                         parinfo[4]["limits"][1]    # 设置恒星大气物理参数第 3 个参数上限 ([Fe/H])
                                         ])

            # 7.14.4.6.2.2 如果 full_output 为 True, 则使用 curve_fit 返回所有信息
            if full_output is True:
                res, res_cov, infodict, errmsg, ier = \
                    curve_fit(min_fun,                                  # 目标函数
                              xdata=[],                                 # x_data 是空列表即可
                              ydata=0.,                                 # y_data 是 min_fun 的理想值
                              p0=value,                                 # 待测参数初始值
                              bounds=bounds,                            # 恒星参数边界值
                              # ftol=np.min([ftol, 1e-4]),              # min_fun 与 ydata 的差异精度
                              ftol=ftol,
                              xtol=xtol,                                # 待测参数的收敛精度
                              # ftol=np.min([ftol, 1e-8]), xtol=1e-10,
                              # ftol=ftol, xtol=1e-10,
                              # options={'maxfev': 10000},
                              full_output=True                          # 打印所有的信息
                              )        
                # 7.14.4.6.2.2.1 计算损失函数
                loss = np.sqrt(np.mean(infodict['fvec'] ** 2))

            # 7.14.4.6.2.3 如果 full_output 为 False, 则使用 curve_fit 返回参数推断均值与标准差
            if full_output is False:
                res, res_cov = \
                    curve_fit(min_fun,                                  # 目标函数
                              xdata=[],                                 # x_data 是空列表即可
                              ydata=0.,                                 # y_data 是 min_fun 的理想值
                              p0=value,                                 # 待测参数初始值
                              bounds=bounds,                            # 恒星参数边界值
                              # ftol=np.min([ftol, 1e-4]),              # min_fun 与 ydata 的差异精度
                              ftol=ftol,
                              xtol=xtol,                                # 待测参数的收敛精度
                              # ftol=np.min([ftol, 1e-8]), xtol=1e-10,
                              # ftol=ftol, xtol=1e-10,              
                              # options={'maxfev': 10000},
                            )
                # 7.14.4.6.2.3.1 计算损失函数
                loss = -9999

            # 7.14.4.6.3 计算标准差, 标准差为协方差对角线矩阵开根号
            result_std = np.diag(res_cov) ** 0.5
            # 7.14.4.6.4 根据误差传播公式计算恒星参数误差
            Rv, Rv_s, Teff, logg, FeH = c0 * (np.exp(signalLog["step"] * res[0]) - 1), c0 * (np.exp(signalLog["step"] * res[1]) - 1),\
                                  np.exp(res[2]), res[3], res[4]
            Rv_err, Rv_s_err, Teff_err, logg_err, FeH_err = result_std[0] * velScale, result_std[1] * velScale, Teff * result_std[2], \
                                                  result_std[3], result_std[4]

            # 7.14.4.6.5 如果 quiet 为 True, 则打印恒星参数推断结果
            # 注意: quiet 为 True 时, full_output 才有效. full_output 为 True 时, plot_fitting 才有效
            # 1) 如果 full_output 为 True, 则打印收敛状态
            # 2) 如果 full_output 为 False, 则不打印收敛状态
            # 3) 如果 plot_fitting 为 True, 则可视化光谱拟合流量残差结果
            if quiet is True:
                print('--------------------------------------------------------------------')
                print('三、推断的恒星大气物理参数')
                print('--------------------------------------------------------------------')
                print("  参数     =   预测值 ±   误差\n"
                        f"1. {'RV':<6}  = {Rv:>8.2f} ± {Rv_err:>6.2f}\n"
                        f"2. {'RV_s':<6}  = {Rv_s:>8.2f} ± {Rv_s_err:>6.2f}\n"
                        f"3. {'Teff':<6}  = {Teff:>8.2f} ± {Teff_err:>6.2f}\n"
                        f"4. {'log g':<6}  = {logg:>8.2f} ± {logg_err:>6.2f}\n"
                        f"5. {'[Fe/H]':<6}  = {FeH:>8.2f} ± {FeH_err:>6.2f}")
                print("--------------------------------------------------------------------")
                if full_output is True:
                    print('--------------------------------------------------------------------')
                    print('四、推断结果收敛状态')
                    print('--------------------------------------------------------------------')
                    print("1. 函数评估次数:", infodict['nfev'])
                    print("2. 最终残差:", infodict['fvec'])
                    # 如果设置 plot_fitting==True, 则可视化光谱拟合结果 (建议针对单个样本检测时使用)
                    if plot_fitting is True:
                        import matplotlib.pyplot as plt
                        from matplotlib.ticker import MaxNLocator
                        plt.style.use('classic')
                        fig = plt.figure(figsize=(30, 8), dpi=300)
                        fig.subplots_adjust(left=0.08, right=0.98, top=0.93, bottom=0.12)
                        plt.plot(np.exp(signalLog["start"] + np.arange(signalLog["data"].shape[0]) * signalLog["step"])[goodPixels],
                                 infodict['fvec'], label="Residual flux")
                        # 如果是 Clean 策略, 则也可视化所 MASK 的异常点
                        if j > 0:
                            # 注意:
                            # 1) 这里的残差不做惩罚, 参考 6.7
                            # 2) 这里的异常点为当前第 j 轮剔除的, 可能与前 j-1 轮的异常点有重复
                            plt.scatter(np.exp(signalLog["start"] + np.setdiff1d(signalLog["goodpix"], goodPixels) * signalLog["step"]),
                                        ((signalLog["data"] - bestfit) / signalLog["err"])[np.setdiff1d(signalLog["goodpix"], goodPixels)],
                                        s=200, c='r', marker='*', label="Outlier")
                        plt.xlabel("Wavelength [Å]", fontsize=30)
                        plt.ylabel("Residual flux", fontsize=30)
                        plt.tick_params(top='on', right='on', which='both', labelsize=25)
                        plt.gca().yaxis.set_major_locator(MaxNLocator(nbins=4))
                        plt.legend(loc="best", fontsize=25)
                        plt.show()
                    # ier=1-4 表示成功, 值越小越好
                    print("3. 拟合状态:", ier)
                    print("4. 拟合结果信息:", errmsg)
                    # 打印状态信息
                    if ier <= 4:
                        print('5. 参数推断成功!')
                    if ier > 4:
                        print('5. 参数推断失败!')
                    print("--------------------------------------------------------------------")

            # 7.14.4.6.6 解析拟合结果
            par_losvd, cmp = uly_fit_pparse(pars=res,              # 待测恒星参数
                                            par_losvd=None,        # 勒让德多项式字典结构
                                            cmp=cmp,               # 字典组件
                                            error=result_std[0:5], # 参数误差
                                            kmoment=kmoment        # 运动学参数数量
                                            )

        # 7.14.4.7 如果待测参数都是已固定的值, 则不需要推断, 直接计算 bestnorm (这里显然不会运行)
        else:
            resc0, bestfit, mpoly, cmp = uly_fitfunc(pars=value,                                              # 待测恒星参数
                                                     kpen=kpen,                                               # 惩罚因子
                                                     cmp=cmp,                                                 # 字典组件
                                                     adegree=adegree,                                         # 勒让德多项式阶数
                                                     signalLog=signalLog,                                     # 光谱字典数据
                                                     mpoly=mpoly,                                             # 勒让德多项式字典结构
                                                     goodpixels=goodPixels,                                   # 有效像素列表
                                                     sampling_function=sampling_function,                     # 插值方法
                                                     voff=voff / velScale,                                    # 速度偏移
                                                     kmoment=kmoment,                                         # 运动学参数数量
                                                     allow_polynomial_reduction=allow_polynomial_reduction    # 是否允许多项式阶数减少
                                                     )
            bestnorm = np.sum(resc0 ** 2)
            if isinstance(parinfo, dict):
                res = parinfo['value']
                error = np.zeros(len(parinfo))
            else:
                error = 0

        # 7.14.4.8 如果达到清理迭代的最大次数 (j 达到 10), 则跳出循环
        if j == nclean:
            break

        # 7.14.4.9 Clean 模式清理异常流量
        # 7.14.4.9.1 保存上一次 Clean 的像素列表, 用于下一次 Clean 的对比
        goodOld = goodPixels.copy()
        # 7.14.4.9.2 计算残差
        resc0, bestfit, mpoly, cmp = uly_fitfunc(pars=value,                                              # 待测恒星参数
                                                 kpen=kpen,                                               # 惩罚因子
                                                 cmp=cmp,                                                 # 字典组件
                                                 adegree=adegree,                                         # 勒让德多项式阶数
                                                 signalLog=signalLog,                                     # 光谱字典数据
                                                 mpoly=mpoly,                                             # 勒让德多项式字典结构
                                                 goodpixels=goodPixels,                                   # 有效像素列表
                                                 sampling_function=sampling_function,                     # 插值方法
                                                 voff=voff / velScale,                                    # 速度偏移
                                                 kmoment=kmoment,                                         # 运动学参数数量
                                                 allow_polynomial_reduction=allow_polynomial_reduction    # 是否允许多项式阶数减少
                                                 )
        # 7.14.4.9.3 计算残差标准差, 并创建残差数组
        rbst0, resc = np.std(resc0, ddof=1), np.zeros(npix)
        # 7.14.4.9.4 获取最佳拟合光谱
        if bestfit is not None:
            bestfit = bestfit
        else:
            bestfit = None
        resc[goodPixels0], bestfit, mpoly, cmp = uly_fitfunc(pars=res,                                                # 待测恒星参数
                                                             kpen=kpen,                                               # 惩罚因子
                                                             cmp=cmp,                                                 # 字典组件
                                                             adegree=adegree,                                         # 勒让德多项式阶数
                                                             signalLog=signalLog,                                     # 光谱字典数据
                                                             mpoly=mpoly,                                             # 勒让德多项式字典结构
                                                             goodpixels=goodPixels,                                   # 有效像素列表
                                                             sampling_function=sampling_function,                     # 插值方法
                                                             voff=voff / velScale,                                    # 速度偏移
                                                             kmoment=kmoment,                                         # 运动学参数数量
                                                             outpixels=goodPixels0,                                   # 输出像素列表
                                                             allow_polynomial_reduction=allow_polynomial_reduction    # 是否允许多项式阶数减少
                                                             )

        """
        第一层裁剪: 剔除明显离群值 (识别宇宙线、错误数据点等明显异常)
        # 注意:
        # 1) 计算未被 MASK 的好像素点的残差 (resc) 的标准差 (rbst_sig)
        # 2) 使用自适应 TGM 模型流量残差阈值 (3 sigma, 4 sigma, 5 sigma, 7 sigma) 识别异常点
        # 3) 判断标准: abs(resc) - modelgrd > clip_level * rbst_sig, 即流量残差能否被模型光谱的波长轻微偏移导致的流量差异解释, 如果不能则该像素点为异常点
        # 4) 如果检测到的点超过 3%, 则提高阈值, 确保不过度剔除
        """
        # 7.14.4.9.5 计算模型中像素的梯度, 以避免由于小像素偏移导致的裁剪问题
        facsh = 0.5 if j == 0 else 0.2
        # 7.14.4.9.6 计算 TGM 模型光谱梯度, 左、右移动 1 个像素
        bestfit_left, bestfit_right= np.roll(bestfit, -1), np.roll(bestfit, 1)
        # 7.14.4.9.6.1 计算左、右光谱梯度绝对值的最大值, facsh 为裁剪因子, noise 为光谱噪声
        modelgrd = np.maximum(np.abs(bestfit - bestfit_left), np.abs(bestfit - bestfit_right)) * facsh / noise
        
        # 7.14.4.9.7 选择大于 3 sigma 或 4,5,7 sigma 的误差 (取决于裁剪分数)
        rbst_sig = np.std(resc[goodPixels], ddof=1)

        # 7.14.4.9.7.1 创建掩码, 选择 "流量残差" 与 "模型流量梯度" 差值大于 3 sigma 的像素
        mask_condition = (np.abs(resc[goodPixels]) - modelgrd[goodPixels]) > (3 * rbst_sig)
        # 7.14.4.9.7.1.1 获取掩码索引
        tmp = np.where(mask_condition)[0].tolist()
        # 7.14.4.9.7.1.2 获取掩码索引数量
        m = len(tmp)
        # 7.14.4.9.7.1.3 设置裁剪级别
        clip_level = 3
        # 7.14.4.9.7.2 根据异常值比例调整裁剪级别, 如果 3 倍标准差剔除了超过 3% 的流量点, 则设置更宽松的剔除级别、即 4 倍标准差
        if m > 0.03 * len(goodPixels):
            # 7.14.4.9.7.2.1 创建掩码
            mask_condition = (np.abs(resc[goodPixels]) - modelgrd[goodPixels]) > (4 * rbst_sig)
            # 7.14.4.9.7.2.2 获取掩码索引
            tmp = np.where(mask_condition)[0].tolist()
            # 7.14.4.9.7.2.3 获取掩码索引数量
            m = len(tmp)
            # 7.14.4.9.7.2.4 设置裁剪级别
            clip_level = 4

            # 7.14.4.9.7.3 根据异常值比例调整裁剪级别, 如果 4 倍标准差剔除了超过 3% 的流量点, 则设置更宽松的剔除级别、即 5 倍标准差
            if m > 0.03 * len(goodPixels):
                # 7.14.4.9.7.3.1 创建掩码
                mask_condition = (np.abs(resc[goodPixels]) - modelgrd[goodPixels]) > (5 * rbst_sig)
                # 7.14.4.9.7.3.2 获取掩码索引
                tmp = np.where(mask_condition)[0].tolist()
                # 7.14.4.9.7.3.3 获取掩码索引数量
                m = len(tmp)
                # 7.14.4.9.7.3.4 设置裁剪级别
                clip_level = 5
                
                # 7.14.4.9.7.4 根据异常值比例调整裁剪级别, 如果 5 倍标准差剔除了超过 3% 的流量点, 则设置更宽松的剔除级别、即 7 倍标准差
                if m > 0.03 * len(goodPixels):
                    # 7.14.4.9.7.4.1 创建掩码
                    mask_condition = (np.abs(resc[goodPixels]) - modelgrd[goodPixels]) > (7 * rbst_sig)
                    # 7.14.4.9.7.4.2 获取掩码索引
                    tmp = np.where(mask_condition)[0].tolist()
                    # 7.14.4.9.7.4.3 获取掩码索引数量
                    m = len(tmp)
                    # 7.14.4.9.7.4.4 设置裁剪级别
                    clip_level = 7

        # 7.14.4.9.8 创建掩码, 并将剔除的 m 个流量点 (clip_level 倍标准差) 的 mask 设置为 0
        # 注意: 从 7.14.4.9.8-7.14.4.9.8.2.4代码有点冗余, 可以优化. 但为了与 IDL 保持一致, 这里保留
        mask = np.zeros(npix)
        # 7.14.4.9.8.1 将 goodPixels0 的 mask 设置为 1
        mask[goodPixels0] = 1
        # 7.14.4.9.8.2 如果 m > 0, 则将剔除的 m 个流量点 (clip_level 倍标准差) 的 mask 设置为 0
        if m > 0:
            # 7.14.4.9.8.2.1 创建最终裁剪级别对应的掩码
            mask_condition = (np.abs(resc[goodPixels0]) - modelgrd[goodPixels0]) > (clip_level * rbst_sig)
            # 7.14.4.9.8.2.2 获取掩码索引
            tmp = np.where(mask_condition)[0].tolist()
            # 7.14.4.9.8.2.3 获取掩码索引
            tmp = np.array(goodPixels0)[tmp]
            # 7.14.4.9.8.2.4 将剔除的 m 个流量点 (clip_level 倍标准差) 的 mask 设置为 0
            mask[tmp] = 0

        """
        第二层裁剪: 根据第一层剔除的坏点, 进一步剔除相邻像素的离群值 (剔除与主要离群点相关的 "边缘" 像素. 处理如天空线残余, 它们通常影响多个相邻像素)
        # 注意:
        # 1) 检查第一层剔除点的直接相邻像素 (左右各一个)
        # 2) 使用较低阈值 (2 sigma) 判断这些相邻点
        # 3) 判断标准: abs(resc[near]) - modelgrd[near] > 2 * rbst_sig
        """
        # 7.14.4.9.9 处理相邻像素
        m2 = 0
        # 7.14.4.9.9.1 如果 m > 0, 则处理相邻像素
        if m > 0:  # 在 2*sig 阈值处清理第一个邻居
            # 7.14.4.9.9.2 获取剔除掉的 m 个流量点的左右邻居
            near = np.concatenate([tmp - 1, tmp + 1])
            # 7.14.4.9.9.3 避免邻居跳出边界
            near = near[(near >= 0) & (near < npix)]
            # 7.14.4.9.9.4 剔除异常邻居, 并更新 tmp
            mask_condition = (np.abs(resc[near]) - modelgrd[near] > 2 * rbst_sig) & (mask[near] == 1)
            # 7.14.4.9.9.5 获取剔除的邻居索引
            nnnn = np.where(mask_condition)[0].tolist()
            # 7.14.4.9.9.6 获取剔除的邻居索引数量
            m2 = len(nnnn)
            # 7.14.4.9.9.7 如果 m2 > 0, 则剔除异常邻居, 并更新 tmp
            if m2 > 0:
                # 7.14.4.9.9.8 将剔除的邻居索引的 mask 设置为 0
                mask[near[nnnn]] = 0
                # 7.14.4.9.9.9 更新 tmp
                tmp = np.concatenate([tmp, near[nnnn]])

        """
        第三层裁剪: 特征完整追踪 (完整追踪和剔除整个光谱连续特征, 如发射线)
        # 注意:
        # 1) 迭代方式追踪光谱特征 (最多 20 次迭代)
        # 2) 每次迭代重新计算 TGM 模型流量残差的标准差 (未被标记为异常)
        # 3) 查找满足五个特定条件的点: 
        #    3.1) 当前未被剔除, mask=1
        #    3.2) 相邻点已被剔除, mask_shifted_back=0
        #    3.3) 当前点与相邻点残差符号相同, ss * resc > 0
        #    3.4) 当前点残差超过标准差, np.abs(resc) - modelgrd > r_sig
        #    3.5) 当前点残差小于相邻已剔除点残差, np.abs(resc) <= np.abs(ss)
        """
        # 7.14.4.9.10 进一步清理
        r_sig = 0
        # 7.14.4.9.10.1 如果 m > 0, 则 mask 异常流量后, 再进一步迭代清理 20 次
        # 注意:
        # 1) mask == 1 表示正常流量
        # 2) mask_shifted_back/mask_shifted_forward 为 mask 向左/右移动 1 个像素后的 mask 表, mask_shifted_back/mask_shifted_forward == 0 表示该像素为正常流量
        # 3) ss 为 resc 向左/右移动 1 个像素后的流量残差 resc, ss * resc > 0 表示该像素的流量为正
        # 4) np.abs(resc) - modelgrd > r_sig 表示 "TGM 模型流量残差" 绝对值与 "TGM 模型流量梯度" 的差值大于 1 倍的 "TGM 模型流量残差" (已剔除异常流量) 的标准差
        # 5) np.abs(resc) <= np.abs(ss) 表示 "TGM 模型流量残差" 绝对值小于等于 "TGM 模型流量残差" 向左/右移动 1 个像素后的绝对值
        if m > 0:
            nclip = 0
            for k in range(1, 21):
                r_sig = np.std(resc[mask == 1], ddof=1)
                # 7.14.4.9.10.2 向前移动
                ss, mask_shifted_back = np.roll(resc, 1), np.roll(mask, 1)
                # 7.14.4.9.10.2.1 满足条件的索引. 
                # 好像素点、但左侧邻居为坏点、且二者同号、且 "TGM 模型流量残差" 与 "流量残差梯度" 之差大于 "TGM 模型流量残差" 标准差、且 "TGM 模型流量残差" 绝对值小于相邻点残差绝对值
                cond1, cond2, cond3, cond4 = (mask == 1) & (mask_shifted_back == 0), (ss * resc > 0), (np.abs(resc) - modelgrd > r_sig), (np.abs(resc) <= np.abs(ss))
                nnn1 = np.where(cond1 & cond2 & cond3 & cond4)[0].tolist()
                mmm1 = len(nnn1)

                # 7.14.4.9.10.3 向后移动
                ss, mask_shifted_forward = np.roll(resc, -1), np.roll(mask, -1)
                # 7.14.4.9.10.3.1 满足条件的索引. 
                # 好像素点、但右侧邻居为坏点、且二者同号、且 "TGM 模型流量残差" 与 "流量残差梯度" 之差大于 "TGM 模型流量残差" 标准差、且 "TGM 模型流量残差" 绝对值小于相邻点残差绝对值
                cond1, cond2, cond3, cond4 = (mask == 1) & (mask_shifted_forward == 0), (ss * resc > 0), (np.abs(resc) - modelgrd > r_sig), (np.abs(resc) <= np.abs(ss))
                nnn2 = np.where(cond1 & cond2 & cond3 & cond4)[0].tolist()
                mmm2 = len(nnn2)

                # 7.14.4.9.10.4 如果 mmm1 > 0, 则将满足条件的索引的 mask 设置为 0
                if mmm1 > 0:
                    mask[nnn1] = 0
                # 7.14.4.9.10.5 如果 mmm2 > 0, 则将满足条件的索引的 mask 设置为 0
                if mmm2 > 0:
                    mask[nnn2] = 0
                # 7.14.4.9.10.6 如果 mmm1 <= 0 且 mmm2 <= 0, 则退出循环
                if (mmm1 <= 0) and (mmm2 <= 0):
                    break
                # 7.14.4.9.10.7 更新 nclip
                nclip += (mmm1 + mmm2)

        # 7.14.4.9.11 更新好像素列表
        if m > 0:
            # 7.14.4.9.11.1 如果 quiet 为 True, 则打印信息.
            # 注意:
            # 1) m + m2 + nclip 为当前第 j 迭代轮剔除的坏像素点, 可能与前 j-1 轮有重复
            if quiet is True:
                print(f'Number of clipped outliers: {m}+{m2}+{nclip}'
                      f' out of {len(goodPixels0)} {rbst_sig}')
            # 7.14.4.9.11.2 更新好像素列表
            goodPixels = np.where(mask == 1)[0].tolist()
        # 7.14.4.9.12 检查是否需要继续迭代. 如果上一次 Clean 的像素列表与当前 Clean 的像素列表相同, 则退出循环
        if np.array_equal(goodOld, goodPixels):
            break

        # 7.14.4.9.13 调整精度
        ftol = 1e-5
        # 7.14.4.9.13.1 如果 j + 2 < nclean, 则设置 ftol = 1e-3
        if j + 2 < nclean:
            ftol = 1e-3
        # 7.14.4.9.13.2 如果 j + 1 < nclean, 则设置 ftol = 1e-4
        elif j + 1 < nclean:
            ftol = 1e-4

        # 7.14.4.9.14 通过统计流量残差的变化, 判断当前优化的效果与剪切操作的影响比例是否合理
        # 根据比例的大小, 决定是从上一个参数位置重新优化, 还是从当前的位置继续优化
        rbst1 = np.std(resc[goodPixels], ddof=1)
        # 7.14.4.9.14.1 如果 rbst_sig ** 2 / (rbst0 * rbst1) < 1.5 且 parinfo 为字典, 则设置 parinfo['value'] = res
        # 注意:
        # 1) rbst_sig 为上一轮 Clean 的 TGM 模型流量与实测光谱流量的残差标准差
        # 2) rbst0 为优化前的 TGM 模型流量与实测光谱流量的残差标准差
        # 3) rbst1 为当前轮 Clean 的 TGM 模型流量与实测光谱流量的残差标准差
        # 4) (rbst_sig ** 2) / (rbst0 * rbst1) 为优化导致的减少与裁剪导致的减少之比
        # if (((rbst_sig ** 2) / (rbst0 * rbst1)) < 1.5) and (isinstance(parinfo, dict)):
        if (((rbst_sig ** 2) / (rbst0 * rbst1)) < 1.5):
            for i in range(5):
                parinfo[i]["value"] = res[i]
        # 7.14.4.9.14.2 优化容差
        elif j + 1 < nclean:
            ftol = 1e-2

    # 7.14.4.9.15 计算推断出 1 条光谱最佳恒星参数的结束时间
    time1 = time.time()
    used_time = time1 - time0
    # 7.14.4.9.15.1 如果 quiet 为 True, 则打印信息
    if quiet is True:
        print("####################################################################")
        print("########################## 耗时：", round(used_time, 2), " s" + " ##########################")
        print("####################################################################")
        print("--------------------------------------------------------------------")

    # 7.14.4.9.16 注意: 参数推断过程中, 我们无法 100% 保证 python 版本与 IDL 版本得到一模一样的结果, 但可得到一致的结果
    # 1) 计算精度差异: IDL 与 Python 计算精度有差异, 特别是在矩阵逆运算时 (如计算勒让德多项式系数).
    # 2) 数值优化方法差异: IDL 与 Python 数值优化方法差异.
    # 3) LASP-CurveFit 会存在恒星参数为 -9999 的异常值 (这些异常值对应的光谱质量较差), 这是由于参数推断失败导致的. 可能包含:
    #    3.1) scipy 的 curve_fit 达到最大迭代次数, 但未收敛.
    #    3.2) 光谱异常, 未开启 Clean 模式. 就算开启 Clean 模式, 也可能没有剔除干净坏流量点, 导致数据处理失败.
    #    3.3) 恒星参数初始值设置不合理, 超出了搜索范围, 导致参数推断失败.
    #    3.4) 矩阵奇异, 逆矩阵运算时导致参数推断失败.
    #    3.5) 对于低质量光谱, 勒让德多项式计算的伪连续谱可能存在负值, 导致参数推断失败.
    #    3.6) 存在未知原因 (我还没有遇到), 导致参数推断失败.
    #    3.7) LASP-Adam-GPU 暂未设置判断推断失败准则.
    # 4) 修正恒星参数误差: 由于 LASP-MPFit 无法对流量不确定性进行无偏估计, 因此得到的参数误差需要改正.
    #    而 LASP-CurveFit 已对流量不确定性进行无偏估计, 并传递到了参数空间, 因此不需要单独修正.
    # 5) 为什么恒星参数误差需要改正: 因为 LASP 计算的卡方函数极小值点, 是基于流量不确定性为 1 的假设.
    #    而实际的流量不确定性与 1 有差异, 这导致恒星参数误差估计值不可靠. 另外, LASP 建立的目标函数, 分母没有考虑流量误差.

    # 7.14.4.9.17 输出 Rv、Teff、log g、[Fe/H]、对应参数误差以及计算时间、流量残差的均方根误差
    return Rv, Teff, logg, FeH, Rv_err, Teff_err, logg_err, FeH_err, used_time, loss
    # 7.14.4.9.18 输出 Rv、Teff、log g、[Fe/H]、对应参数误差以及计算时间、流量残差的均方根误差、Clean 次数, 剔除的异常点数量
    # return Rv, Teff, logg, FeH, Rv_err, Teff_err, logg_err, FeH_err, used_time, loss, j, len(np.setdiff1d(signalLog["goodpix"], goodPixels))
