# -*- coding: utf-8 -*-
# @Time    : 11/12/2024 15.18
# @Author  : ljc
# @FileName: uly_fit_lin.py
# @Software: PyCharm


# 1. 简介
"""
 Python conversion of the IDL uly_fit_lin.pro .
目的：
    对 TGM 光谱进行采样插值、降低 TGM 光谱分辨率、乘以勒让德多项式, 用于与待测光谱匹配.
函数：
    1) uly_fit_fractsolve
    2) uly_fit_lin_mulpol
    3) uly_fit_lin_weight
    4) uly_fit_lin
解释：
    1) uly_fit_fractsolve 函数: 返回 TGM 光谱与待测光谱的权重, 由 uly_fit_lin_weight 函数调用;
    2) uly_fit_lin_mulpol 函数: 建立勒让德乘法多项式, 返回 TGM 光谱 (乘以多项式)、更新后的 mpoly;
    3) uly_fit_lin_weight 函数: 返回乘以勒让德多项式 * 权重后的 TGM 光谱、更新后的 cmp、None 作为可加性多项式;
    4) uly_fit_lin 函数: 返回乘以勒让德多项式、加权、降低分辨率后的 TGM 光谱、勒让德多项式.
"""


# 2. 调库
import numpy as np
from legendre_polynomial.mregress import mregress
from uly_tgm_eval.uly_tgm_eval import uly_tgm_eval
from resolution_reduction.convol import convol
import warnings
warnings.filterwarnings("ignore")


# 3. 计算模型中每个分量 cmp 的权重系数 (LAMOST 的 cmp 是 1 个字典, 因此 cmp 的长度记为 1)
def uly_fit_fractsolve(a=None, b=None, npoly=None, bounds=None) -> np.ndarray:

    """
        计算模型中分量 cmp 的权重系数.
        注意: 
            1) LASP-CurveFit 中读取 LASMOT 光谱的函数 uly_spect_read_lms 支持流量中值归一化, 因此该函数意义不大, 但为了与 LASP-MPFit 一致, 这里暂留.
            2) 该函数的目标是: 最小化 ||b - w*a||², 使得模型光谱与实测光谱量纲尽可能一致, 解析解: w = (a^T * b) / (a^T * a).

        输入参数:
        -----------
        a:
          1 个向量, 直接用最简单的线性拟合公式计算权重.
          注意: 如果有另外需求, 请参考 ulyss 原 IDL 代码.
        b:  
          待测光谱.
        npoly: 无用参数, 但为了与 LASP-MPFit 一致, 这里暂留，
              多项式基的个数.
        bounds: 无用参数, 但为了与 LASP-MPFit 一致, 这里暂留.
               权重系数的上下限.

        输出参数:
        -----------
        soluz:
              该函数返回一个权重数组, 它被 uly_fit_lin_weight 调用.
    """

    # 3.1 将 a 和 b 转换为 numpy 数组, 并重塑为列向量
    a = np.array(a).reshape(-1, 1)
    b = np.array(b).reshape(-1, 1)

    # 3.2 计算权重
    soluz = np.array([np.sum(a * b) / np.sum(a ** 2)])

    # 3.3 返回权重
    return soluz


# 4. 使用勒让德多项式构建流量改正因子, 改正模型光谱与实测光谱形状差异
def uly_fit_lin_mulpol(bestfit=None, cmp=None, mpoly=None, SignalLog=None, goodPixels=None, polpen=None, allow_polynomial_reduction=False) -> tuple[np.ndarray, dict]:

    """
        使用勒让德多项式构建伪谱 (流量改正因子), 改正模型光谱与实测光谱形状差异.

        输入参数:
        -----------
        bestfit:
                TGM 插值光谱 * 伪光谱 (光谱改正因子) ≈ 实测光谱.
        cmp:
            为字典格式. 存储 TGM、待测参数初始值、勒让德多项式默认值等信息;
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
        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).
        SignalLog:
                  包含待测光谱的结构体 (见 CMP), 模型光谱和待测光谱的波长必须是 ln 对数化的, 相关标签在 uly_spect_alloc 文档中有描述.
        goodpixels:
                   好像素点的索引值.
        polpen:
               乘法多项式的偏置水平. 此关键词可用于减少乘法多项式中不重要项的影响.
              注意: 
              1) 默认情况下不应用偏置. 如果某些系数的绝对值小于 polpen 倍的统计
                 误差, 这些系数会通过因子 (abs(coef)/(polpen*err))^2 被抑制.
              2) 该功能仅在 mdegree>0 时有效, polpen=2 是一个合理的选择.
        allow_polynomial_reduction:
                                   是否允许多项式阶数减少. 默认值为 False, 即不允许多项式阶数减少.
                                   如果设置为 True, 则允许多项式阶数减少.
                                   注意: 伪光谱 (流量改正因子) 不应该小于 0, 如果小于 0, 提供两种处理方法:
                                   1) 为避免伪光谱谱存在负值, 循环减少多项式阶数, LASP-MPFit 版本采用这种方法.
                                   2) 直接认为该光谱质量较差, 参数推断失败, LASP-CurveFit 默认采用这种方法 (因为伪谱为负值, 大概率由于光谱流量存在负值).
                                   3) 当前 LASP-Adam-GPU 版本不提供该功能, 如果 LASP-Adam-GPU 流量改正因子为负值, 同样得到最佳参数 (不输出 -9999), 因此这类数据所测参数可能不够可靠.

        输出参数:
        -----------
        bestfitp:
                 TGM 光谱 (乘以伪谱, 即乘以流量改正因子后的模型光谱).
        mpoly:
              更新后的 mpoly 字典信息.
    """

    # 4.1 获取待测光谱的好像素点流量、好像素点流量误差
    # 注意: 
    # 1) LASP 默认好像素点流量误差为 1 (没有使用倒方差), 因此使用 MPFit 推断的参数误差需改正 (需无偏估计计算光谱流量的误差)
    galaxy, noise = SignalLog["data"][goodPixels], SignalLog["err"][goodPixels]
    # 4.1.1 如果光谱流量、或好像素点流量误差全为非正值, 则认为该光谱质量较差, 参数推断失败
    # 流量没有正数, 则抛出异常
    if (np.max(galaxy) <= 0) or (np.max(noise) <= 0):
        raise ValueError("The spectrum flux or error has all negative value, please check the spectrum quality!")

    # 4.2 计算 TGM 插值光谱
    # 注意: 
    # 1) 伪谱 (流量改正因子): (1, c1, c2, ..., cn) * (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)) = (1, c1 * P1(lambda_i), c2 * P2(lambda_i), ..., cn * Pn(lambda_i))
    # 2) 伪谱的初始值为单位向量, 即 (1, 1, 1, ..., 1)
    # 2) 待解值是勒让德多项式的系数 (1, c1, c2, ..., cn), 需使用 mregress 函数求解
    # 3) 实测光谱 ≈ 伪谱 (流量改正因子) * TGM 模型光谱 ≈ (1, c1, c2, ..., cn) * (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)) * TGM 模型光谱
    # 4.2.1 bestfitp 为 TGM 模型光谱, bestfit 为乘以了伪谱的 TGM 模型光谱, mpoly["poly"] 的形状是 (流量维度, 1)
    bestfitp = bestfit / mpoly["poly"].reshape(-1, 1)

    # 4.2.2 cmul = (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)) * TGM 光谱, mpoly["leg_array"] 的形状是 (流量维度, 多项式项数)
    cmul = mpoly["leg_array"] * bestfitp
    # 4.2.3 执行多元线性回归拟合, 解出回归系数 coefs_pol, 即 (1, c1, c2, ..., cn) 
    # 注意: mregress 为自定义函数, 与 LASP-MPFit 一致, coefs_pol 的形状是 (多项式项数,)
    coefs_pol, inv = mregress(cmul[goodPixels, :],   # 伪谱 (流量改正因子) * TGM 光谱
                              galaxy,                # 待测光谱好像素点流量
                              measure_errors=noise   # 待测光谱好像素点流量误差
                              )
    # 4.2.3.1 如果 polpen 不为 None, 则使用偏置水平
    # 注意:
    # 1) 如果 polpen 为 None, 则不使用偏置水平, LASP-MPFit 默认不使用偏置水平, 因此代码块 4.2.3.1 可以删除
    # 2) 我们保留了与 LASP-MPFit 一致的代码, 详情可参考 IDL 代码
    # if polpen is not None:
    #     penal_pol = np.abs(coefs_pol) / (polpen * inv["sigma"])
    #     indices  = np.where(penal_pol[1:] < 1)[0]
    #     if len(indices) > 0:
    #         pen = 1 + indices
    #         coefs_pol[pen] *= penal_pol[pen] ** 2 

    # 4.2.4 更新伪谱 (流量改正因子) mpoly["poly"]
    # 注意: IDL 与 Python 索引取值有差异 (Python 不包括最后一个索引, 但 IDL 包括最后一个索引)
    # 4.2.4.1 伪谱 (流量改正因子): (1, c1, c2, ..., cn) * (1, P1(lambda_i), P2(lambda_i), ..., Pn(lambda_i)) = (1, c1 * P1(lambda_i), c2 * P2(lambda_i), ..., cn * Pn(lambda_i))
    mpoly["poly"] = mpoly["leg_array"][:, 0: mpoly["lmdegree"] + 1]  @ np.transpose(coefs_pol).reshape(-1, 1)
    minmaxpol = [np.min(mpoly["poly"][goodPixels]), np.max(mpoly["poly"][goodPixels])]

    # 4.3 伪谱 (流量改正因子) 不应该小于 0, 如果小于 0, 提供两种处理方法:
    # 1) 为避免伪谱 (流量改正因子) 存在负值, 循环减少多项式阶数, LASP-MPFit 版本采用这种方法
    # 2) 直接认为该光谱质量较差, 参数推断失败, LASP-CurveFit 版本默认采用这种方法 (因为伪谱为负值, 大概率由于光谱流量存在负值)
    # 4.3.1 如果伪谱 (流量改正因子) 最大值小于 0, 则认为该光谱质量较差, 参数推断失败.
    if minmaxpol[1] <= 0:
        raise ValueError("Pseudo-continuum all flux values are less than 0, please check the spectrum quality!")
    # 4.3.2 方法1) allow_polynomial_reduction 为 True 时, 如果伪谱 (流量改正因子) 最小值小于 0, 则减少多项式项数
    if (minmaxpol[0] < 0) & (minmaxpol[1] > 0) & (allow_polynomial_reduction):
        # 4.3.2.1 获取勒让德多项式最大阶数
        lmd = mpoly["lmdegree"]
        # 4.3.2.2 如果伪谱小于 0, 则减少多项式项数
        while (np.min(mpoly["poly"]) <= 0) & (lmd > 0):
            lmd = lmd - 1
            # 4.3.2.3 执行多元线性回归拟合, 计算回归系数 coefp
            coefp, inv = mregress(cmul[goodPixels, 0: lmd + 1],  # 伪谱 (流量改正因子) * TGM 光谱
                                  galaxy,                        # 待测光谱好像素点流量
                                  measure_errors=noise           # 待测光谱好像素点流量误差
                                 )
             # 4.3.2.3.1 如果 polpen 不为 None, 则使用偏置水平
             # 注意:
             # 1) 如果 polpen 为 None, 则不使用偏置水平, LASP 默认不使用偏置水平, 因此代码块 4.3.2.3.1 可以删除
             # 2) IDL 中有处代码错误, coefs_pol[pen] *= penal_pol[pen]^2 应改为 coefp[pen] *= penal_pol[pen]**2
             # 3) 除此之外, 我们保留了与 LASP-MPFit 一致的代码, 详情可参考 IDL 代码
            if polpen is not None:
                penal_pol = np.abs(coefp) / (polpen * inv["sigma"])
                indices  = np.where(penal_pol[1:] < 1)[0]
                if len(indices) > 0:
                    pen = 1 + indices
                    coefp[pen] *= penal_pol[pen] ** 2 
            # 4.3.2.4 更新伪谱 (流量改正因子) mpoly["poly"], 阶数为 lmd
            mpoly["poly"] = mpoly["leg_array"][:, 0: lmd + 1]  @ np.transpose(coefp).reshape(-1, 1)
        # 4.3.2.5 更新回归系数 coefs_pol
        # coefs_pol = np.zeros_like(coefs_pol) + coefp
        coefs_pol = coefp
    # 4.3.3 方法2) allow_polynomial_reduction 为 False 时, 如果伪谱 (流量改正因子) 最小值小于 0, 也认为该光谱质量较差, 参数推断失败
    if (minmaxpol[0] < 0) & (minmaxpol[1] > 0) & (not allow_polynomial_reduction):
        raise ValueError("Pseudo-continuum is less than 0, please check the spectrum quality!")

    # 4.4 乘以伪谱 (流量改正因子) 的 TGM 光谱 bestfitp
    bestfitp = bestfitp * mpoly["poly"]

    # 4.5 更新勒让德多项式系数 mpoly["mpolcoefs"]
    mpoly["mpolcoefs"] = coefs_pol
    
    # 4.6 返回乘以伪谱 (流量改正因子) 的 TGM 光谱 (bestfitp), 更新后的勒让德多项式系数 mpoly["mpolcoefs"]
    return bestfitp, mpoly


# 5. 计算流量比
def uly_fit_lin_weight(SignalLog=None, goodpixels=None, tmp=None, mpoly=None, adegree=None, lum_weight=None, cmp=None, addcont=None) -> tuple[np.ndarray, dict, np.ndarray]:

    """
        确定每个成分的权重, LASP 只使用了 1 个成分 (cmp).

        输入参数:
        -----------
        SignalLog:
                  包含待测光谱的结构体 (见 CMP), 模型光谱和待测光谱的波长必须是 ln 对数化的, 相关标签在 uly_spect_alloc 文档中有描述.
        goodpixels:
                   好像素点的索引值.
        tmp:
            TGM 模板.
        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).
        adegree:
                加性多项式的阶数, LASP 设置为 -1 (不使用加性多项式), 我认为该函数的功能可以消除理论光谱与实测光谱的差异, 后续可用于理论光谱模拟器推断实测光谱参数.
        lum_weight:
                   计算流量权重及其误差. 如果设置, 计算 e_weight 和 l_weight.
        cmp:
            TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
        addcont:
                加性多项式数组.

        输出参数:
        -----------
        bestfit:
                乘以流量改正因子、权重后的 TGM 光谱.
        cmp:
            结果以 cmp["weight"] 返回. 如果设置 lum_weight, 则包含 e_weight, l_weight.
        addcont:
                加性谱 (数组).
    """

    # 5.1 获取成分 (组件) 的数量 (LASP 设置为 1)、待拟合的流量点数量. 组件是什么?????? 我没有理解透彻
    npix, n_cmp = SignalLog["data"].shape[0], 1

    # 5.2 获取待拟合的好像素点的流量、待拟合的好像素点的流量误差, 需注意: LASP 没有使用 LAMOST 流量的真实误差
    galaxy, noise = np.double(SignalLog["data"][goodpixels]), np.double(SignalLog["err"][goodpixels])

    # 5.3 实测光谱 ≈ 伪谱 (流量改正因子) * TGM 模型光谱
    c = np.double(mpoly["poly"]).reshape(-1, 1) * tmp
    """
    # 5.3.1 如果有加性多项式, 在前面添加 addpol*mulpol.
    # 由于 LASP-MPFit 中 adegree=-1. 如有需要, 请参考 IDL 代码.
    # if adegree >= 0:
    #     c1 = np.double(mpoly["poly"]).reshape(-1, 1)
    #     x = 2.0 * np.arange(npix, dtype=np.float64) / npix - 1.0
    #     for j in range(adegree + 1):
    #         c1[:, j] *= legendre(x, j)
    #     c = np.hstack([c1, c])
    """

    # 5.4 好像素点流量
    a = c[goodpixels, :]
    for j in range(0, (adegree + 1) + n_cmp - 1 + 1):
        # 5.4.1 好像素点流量除以好像素点流量误差, LAMOST 设置的 noise 为 1, 因此不使用倒方差
        a[:, j] /= noise
    """
    # 5.4.2 Tikhonov 正则化部分, LASP-MPFit 没有使用, 因此不做调试. 如有需要, 请参考原始 IDL 代码.
    lambda_reg = 0  # 正则化权重为 0 表示不使用正则化
    # if lambda_reg > 0:
    #     ng = len(goodpixels)
    #     nrel = ng + n_cmp - 1
    #     aa = np.zeros((nrel, n_cmp))
    #     aa[0, 0] = a[:, 0]
    #     a = aa
    #     cr = np.zeros(c.shape[1])
    #     cr[adegree + 1:adegree + 3] = lambda_reg * np.array([-1, 1])  # 1st order regularization
    #     for k in range(ng, nrel):
    #         a[k, :] = np.roll(cr, k - ng)
    #     goodpixels = np.concatenate([goodpixels, np.arange(n_cmp - 1) + npix])
    #     galaxy = np.concatenate([galaxy, np.zeros(n_cmp - 1)])
    #     noise = np.concatenate([noise, np.ones(n_cmp - 1)])
    """

    # 5.5 拟合中可加多项式的阶数, LASP 设置的 adegree=-1, 所以 npoly=0
    npoly = (adegree + 1)

    # 5.6 处理拟合
    # 5.6.1 如果 adegree=-1, 则计算实测光谱与 TGM 光谱的权重 (倍数关系)
    if adegree == -1:
        # 5.6.1.1 计算权重, 计算实测光谱与 TGM 光谱的权重 (倍数关系, 这一步可以不要, 影响不大)
        cmp["weight"] = uly_fit_fractsolve(a,                     # 模型光谱好像素点流量, 已乘以伪谱 (流量改正因子) 的 TGM 光谱
                                           galaxy / noise,        # 待测光谱好像素点流量/流量误差
                                           npoly,                 # 可加多项式阶数, LASP 设置为 0
                                           bounds=cmp["lim_weig"] # 权重限制
                                           )
        # 5.6.1.2 乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱
        # 注意: "权重" 进一步改正了模型光谱与实测光谱的倍数差异
        bestfit = c * cmp["weight"]
        # 5.6.1.3 LASP-MPFit 设置的 adegree=-1, 因此可加性多项式为 0
        addcont = np.zeros_like(bestfit)
    # 5.6.2 如果设置加性多项式, 则计算实测光谱与 TGM 光谱的权重
    else:
        # 5.6.2.1 计算权重, 计算实测光谱与 TGM 光谱的权重 (倍数关系)
        wght = uly_fit_fractsolve(a,                       # 模型光谱好像素点流量, 已乘以伪谱 (流量改正因子) 的 TGM 光谱
                                  galaxy / noise,          # 待测光谱好像素点流量/流量误差
                                  npoly,                   # 可加多项式的阶数
                                  bounds=cmp["lim_weig"]   # 权重限制
                                  )
        # 5.6.2.2 计算乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱
        # bestfit = c @ wght
        bestfit = c * wght
        # 5.6.2.3 更新 cmp["weight"]
        cmp["weight"] = wght[adegree + 1:]
        addcont = c[:, 0: adegree + 1] @ wght[0: adegree + 1]

    # 5.7 检查权重无效值, 并把无效值设置为 0 (注意: cmp["weight"] 为 1 个比例系数)
    # 5.7.1 检查是否无穷大无,
    # 注意:
    # 1) 我仅根据 IDL 更改, Python 可能会有更高效的方法, 留给后人更新吧, 比如: 
    # if (np.isfinite(cmp["weight"])) | (np.abs(cmp["weight"]) > 1e36) | (np.abs(cmp["weight"]) < 1e-36):
        # cmp["weight"] = 0
    nan_idx = np.where(~np.isfinite(cmp["weight"]) | (np.abs(cmp["weight"]) > 1e36))[0].tolist()
    if len(nan_idx) > 0:
        cmp["weight"] = 0
    # 5.7.2 处理接近零的值
    zeros_idx = np.where(np.abs(cmp["weight"]) < 1e-36)[0].tolist()
    if len(zeros_idx) > 0:
        cmp["weight"] = 0

    # 5.8 更新 cmp["e_weight"] 和 cmp["l_weight"]
    # 5.8.1 如果 lum_weight 为 None, 并且 cmp["weight"] 不为空
    if (lum_weight is None) and (cmp["weight"] != ""):
        # 5.8.1.1 获取 cmp["weight"] 不为 0 的索引
        w_notnull = np.where(cmp["weight"] != 0)[0].tolist()
        # 5.8.1.2 初始化 cmp["e_weight"] 和 cmp["l_weight"]
        cmp["e_weight"], cmp["l_weight"] = np.zeros(1), np.zeros(1)
        # 5.8.1.3 如果 w_notnull 不为空
        if len(w_notnull) != 0:
            inv = 0
            if adegree > -1:
                w_nn = np.concatenate([np.arange(adegree + 1), w_notnull + adegree + 1])
            else:
                w_nn = w_notnull
            ww, inv = mregress(c[goodpixels][:, w_nn],  # 模型光谱好像素点流量, 已乘以伪谱 (流量改正因子) 的 TGM 光谱
                               galaxy,                  # 待测光谱好像素点流量
                               measure_errors=noise,    # 待测光谱好像素点流量误差
                               inv=inv                  # 接收中间结果
                               )
            cmp["e_weight"], cmp["l_weight"] = inv["sigma"][adegree + 1:], (np.sum(tmp, axis=0) * cmp["weight"]) / npix
    else:
        # 5.8.2 如果 lum_weight 不为 None, 并且 cmp["weight"] 为空
        cmp["e_weight"], cmp["l_weight"] = 0, 1

    # 5.9 返回乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱, 更新后的 cmp、None 作为可加性多项式
    # 5.9.1 如果 adegree == -1, 则返回乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱, 更新后的 cmp、None 作为可加性多项式
    if adegree == -1:
        return bestfit, cmp, None
    # 5.9.2 如果 adegree != -1, 则返回乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱, 更新后的 cmp、可加性多项式
    else:
        return bestfit, cmp, addcont


# 6. 拟合线性系数, 并返回最佳拟合结果
def uly_fit_lin(adegree=None, par_losvd=None, cmp=None, SignalLog=None, goodPixels=None, voff=None, mpoly=None, sampling_function=None,
                addcont=None, lum_weight=None, polpen=None, cache=None, modecvg=None, allow_polynomial_reduction=False) -> tuple[np.ndarray, dict]:

    """
        该函数调整乘法和加法连续项 (LASP 只使用乘法连续项), 使用最小二乘法计算最佳拟合 (勒让德多项式系数), 并返回最佳拟合结果.

        输入参数:
        -----------
        adegree:
                用于修正 TGM 模板光谱形状的加法勒让德多项式的阶数. 在拟合过程中, 默认不使用任何加法多项式.
                如果要禁用加法多项式, 请设置 adegree=-1, LASP 设置为 -1.
        par_losvd:
                  losvd 参数, 数组中包含 0 到 6 个数, 顺序为: [cz, sigma, h3, h4, h5, h6], 其中
                  cz 和 sigma 为像素值. LASP 仅使用 cz、sigma 降低 TGM 光谱分辨率, 并计算 Rv.
        cmp:
            TGM 字典, 存储 TGM、待测参数初始值、勒让德多项式默认值等信息.
            cmp 是 1 个字典而不是字典列表, 因此 cmp 的长度记为 1.
        SignalLog:
                  包含待测光谱的结构体 (见 CMP), 模型光谱和待测光谱的波长必须是 ln 对数化的, 相关标签在 uly_spect_alloc 文档中有描述.
        goodpixels:
                   好像素点的索引值.
        voff:
             待分析光谱与模型之间的速度偏移, 单位为 km/s.
        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).
        sampling_function:
                         插值方法. 可输入 "splinf", "cubic", "slinear", "quadratic", "linear". 默认使用 "linear" 插值方法.
        addcont:
                加法多项式数组.
                注意:
                1) 在拟合过程中, 默认不使用任何加法多项式.
        lum_weight:
                   计算流量权重及其误差. 如果设置, 计算 e_weight 和 l_weight.
        polpen:
               乘法多项式的偏置水平. 此关键词可用于减少乘法多项式中不重要项的影响.
               默认情况下不应用偏置. 如果某些系数的绝对值小于 polpen 倍的统计
               误差, 这些系数会通过因子 (abs(coef)/(polpen*err))^2 被抑制.
               该功能仅在 mdegree>0 时有效, polpen=2 是一个合理的选择.
        cache:
              初始化缓存.
        modecvg:
                指定拟合方法的收敛模式. LASP-MPFit 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码.
                注意:
                1) modecvg=0 (默认选项), 这是最快的, 但如果解错失, 问题可能发生.
                2) modecvg=1 (每次迭代仅完成一次), 为每个 LM 迭代计算导数.
                3) modecvg=2, uly_fit_lin 始终收敛, 但速度较慢.
        allow_polynomial_reduction:
                                   是否允许多项式阶数减少. 默认值为 False, 即不允许多项式阶数减少.
                                   如果设置为 True, 则允许多项式阶数减少.
                                   注意: 伪光谱 (流量改正因子) 不应该小于 0, 如果小于 0, 提供两种处理方法:
                                   1) 为避免伪光谱谱存在负值, 循环减少多项式阶数, LASP-MPFit 版本采用这种方法.
                                   2) 直接认为该光谱质量较差, 参数推断失败, LASP-CurveFit 默认采用这种方法 (因为伪谱为负值, 大概率由于光谱流量存在负值).
                                   3) 当前 LASP-Adam-GPU 版本不提供该功能, 如果 LASP-Adam-GPU 流量改正因子为负值, 同样得到最佳参数 (不输出 -9999), 因此这类数据所测参数可能不够可靠.

        输出参数:
        -----------
        bestfit:
                乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱.
        cmp:
            更新后的 cmp 字典.
    """

    # 6.1 获取待测光谱流量维度、成分数量 (也就是 cmp 的数量)、好像素点流量、好像素点流量误差
    npix, n_cmp, galaxy, noise = SignalLog["data"].shape[0], 1, SignalLog["data"][goodPixels], SignalLog["err"][goodPixels]

    # 6.2 默认勒让德多项式初始字典结构 (lmdegree: 勒让德多项式的最大阶数、mpolcoefs: 勒让德多项式系数、 poly: 伪谱, 即流量改正因子)
    if mpoly is None:
        mpoly = {'lmdegree': 0, 'mpolcoefs': np.zeros(1), 'poly': np.ones(npix)}

    # 6.3 检查待测光谱流量噪声 (误差) 是否严格为正
    if np.min(noise) <= 0:
        raise ValueError('NOISE vector must be strictly positive!')

    # 6.4 如果没有提供, 初始化缓存 (para: 3 个大气参数、 spectra: 初始模型光谱流量)
    if cache is None:
        cache = {'para': [-1e35, -1e35, -1e35] * n_cmp, 'spectra': np.zeros((npix, n_cmp))}

    # 6.5 初始化 TGM 模型谱
    models = np.zeros((npix, n_cmp))

    # 6.6 评估所有组件
    # 1) cmp 为 1 维, 即 1 个字典
    # 2) cmp['para'] 包含 Teff、log g、[Fe/H] 3 个参数, 因此 len(cmp['para'])=3
    cached = 0
    if len(cmp['para']) > 0:
        # 6.6.1 如果大气物理参数值达到极限小, 则 TGM 光谱流量为 0
        para_values = np.array([param['value'] for param in cmp['para']])
        if np.all(para_values == cache['para']):
            models[:, 0] = cache['spectra'][:, 0]
            cached = 1

    # 6.7 如果大气物理参数值正常, 则调用 TGM 模型, 生成 TGM 光谱
    # 注意: LASP 使用的是 uly_tgm_eval, 前两个(SSP, STAR), 如果需要请参考原始 IDL 代码
    if cached == 0:
        # 6.7.1 如果 cmp['eval_fun'] == 'SSP', 则调用 uly_ssp_interp 函数, 生成 TGM 光谱
        if cmp['eval_fun'] == 'SSP':
            # 6.7.1.1 调用 uly_ssp_interp 函数, 生成 TGM 光谱
            mdl = uly_ssp_interp(cmp['eval_data'], cmp['para']['value'])
        # 6.7.2 如果 cmp['eval_fun'] == 'STAR', 则 mdl = cmp['eval_data']
        elif cmp['eval_fun'] == 'STAR':
            mdl = cmp['eval_data']
        else:
            # 6.7.3 指定 para_values 时, 调用 uly_tgm_eval 函数生成 TGM 光谱
            para_values = np.array([param['value'] for param in cmp['para']])
            # 6.7.3.1 调用 uly_tgm_eval 函数生成 TGM 光谱
            mdl = uly_tgm_eval(eval_data=cmp['eval_data'],         # 模型光谱系数
                               para=para_values,                   # 大气物理参数
                               sampling_function=sampling_function # 插值方法
                               )

        # 6.7.4 更新 TGM 模型光谱 (感觉在 Python 中无用, 但为了与 IDL 一致, 这里进行定义)
        # 注意: IDL 索引取值包括最后一个、Python 索引取值不包含最后一个, 这是 Python 与 IDL 的差异
        nmd = np.min([npix, len(mdl)])
        models[0: nmd, 0] = mdl[0: nmd]
        para_values = np.array([param['value'] for param in cmp['para']])

        # 6.7.5 更新 cache 中的 Teff、log g、[Fe/H] 值, 这段代码 (6.7.5) 块其实可以删除!
        # len(cmp['para']) 为 3
        if len(cmp['para']) > 0:
            # 6.7.5.1 更新 cache 中的 Teff、log g、[Fe/H] 值
            cache['para'] = para_values
        # 6.7.5.2 更新 cache 中的 模型光谱值
        cache['spectra'][:, 0] = models[:, 0]

    # 6.8 使用 losvd 参数对 TGM 光谱进行卷积、降低分辨率到 LAMOST 光谱
    # 注意: 
    # 1) LASP-MPFit 使用了 losvd 的前两个参数
    # 2) scipy 优化器迭代优化 losvd 参数. 因此每一轮迭代, 降低分辨率的参数不一样, 最终使用与 LAMOST 光谱卡方值最小的参数
    if par_losvd is not None:
        # 6.8.1 高斯核函数的均值
        maxvel = np.abs(par_losvd[0])
        # 6.8.2 高斯核函数的标准差
        maxsig = par_losvd[1]
        # 6.8.3 使用 5 倍标准差范围的长度
        dx = np.ceil(np.abs(voff) + np.abs(maxvel) + 5.0 * maxsig)
        # 6.8.4 高斯核的半宽, 不超过光谱长度的一半
        dx = np.min([dx, (npix - 1.0) / 2.0])
        # 6.8.5 高斯核函数的总长度, 使用 "奇数" 长度
        n = int(2 * dx + 1)
        # 6.8.6 x 关于中心位置互为相反数 (对称)
        x = dx - np.arange(n)
        # 6.8.7 保持与 LASP-MPFit 类似的初始
        losvd = np.empty(n, dtype=np.float64)
        # 6.8.8 速度加上偏移
        vel = voff + par_losvd[0]
        # 6.8.9 高斯核的标准差 (为啥标准差不小于 0.1)
        sigma_pix = np.max([par_losvd[1], 0.1])
        # 6.8.10 计算高斯核的值 exp((-(x-mu)/sigma) ** 2)
        w = (x - vel) / sigma_pix
        w2 = w * w
        # 6.8.11 用 0 替换大于 5*sigma 的值, 然后计算高斯核的值 exp((-(x-mu)/sigma) ** 2)
        wlarge = np.where(np.abs(w) > 5.0)[0].tolist()
        wnorm = np.where(np.abs(w) <= 5.0)[0].tolist()
        if len(wlarge) > 0:
            losvd[wlarge] = 0.
        # 6.8.12 赋值 losvd, 即核函数的值, 用于降低分辨率
        losvd[wnorm] = np.exp(-0.5 * w2[wnorm]) / (np.sqrt(2.0 * np.pi) * sigma_pix)

        # 6.8.13 处理厄米特多项式, LASP 没有使用 h3, ... , h6
        # 6.8.13.1 厄米特多项式的个数
        # 注意: LASP 没有使用 h3, ... , h6, 因此 nherm=0
        nherm = len(par_losvd) - 2
        # 6.8.13.2 如果 nherm > 0, 则计算厄米特多项式
        if nherm > 0:
            # 6.8.13.2.1 如果 nherm=1, 则计算厄米特多项式 h3
            poly = 1.0 + par_losvd[2] / np.sqrt(3.0) * (w * (2.0 * w2 - 3.0))
            # 6.8.13.2.2 如果 nherm > 1, 则计算厄米特多项式  h4
            if nherm > 1:
                poly += par_losvd[3] / np.sqrt(24.0) * (w2 * (4.0 * w2 - 12.0) + 3.0)
            # 6.8.13.2.3 如果 nherm > 2, 则计算厄米特多项式 h5
            if nherm > 2:
                poly += par_losvd[4] / np.sqrt(60.0) * (w * (w2 * (4.0 * w2 - 20.0) + 15.0))
            # 6.8.13.2.4 如果 nherm > 3, 则计算厄米特多项式 h6
            if nherm > 3:
                poly += par_losvd[5] / np.sqrt(720.0) * (w2 * (w2 * (8.0 * w2 - 60.0) + 90.0) - 15.0)
            # 6.8.13.2.5 乘以厄米特多项式
            losvd *= poly

        # 6.8.14 归一化核函数   
        losvd /= np.sum(losvd)

        # 6.8.15 使用 losvd 参数卷积 TGM 光谱, 降低 TGM 光谱分辨率到 LAMOST
        # 注意: 使用自定义 convol 函数保证与 LASP-MPFit 结果一致
        tmp = convol(models,   # 高分辨率 TGM 模型光谱
                     losvd     # 核函数, 用于降低 TGM 光谱分辨率到 LAMOST
                     )
    else:
        tmp = models.copy()
    

    # 6.9 迭代循环设置
    """
    注意: 对于 LASP 来说, 只使用了 1 个成分 (cmp), 因此原始 IDL 代码有点冗余. 这里提供两种方案:
    1) 方案 1 (6.9.1): 删除冗余代码, 参数推断结果几乎不影响, 如 Teff 可能受数值精度影响产生 0.1 K 左右差异
    2) 方案 2 (6.9.2): 保留与 IDL 一致的代码, 目前也使用该方案. 详细见 IDL 代码
    3) 无论哪种方案, 均先计算流量比消除模型与待测光谱的比例差异, 然后再乘以伪谱 (流量改正因子) 消除模型光谱与待测光谱的形状差异
    4) 如何取舍, 取决于个人喜好
    """
    # 6.9.1 方案 1
    # bestfitw, cmp, addcont = uly_fit_lin_weight(SignalLog,              # 待测光谱字典结构
    #                                             goodPixels,             # 有效像素列表
    #                                             tmp,                    # 降低分辨率后的 TGM 模型光谱
    #                                             mpoly,                  # 勒让德多项式字典结构
    #                                             adegree=adegree,        # 加型多项式阶数
    #                                             lum_weight=lum_weight,  # 计算流量权重及其误差
    #                                             cmp=cmp,                # TGM 模型字典结构
    #                                             addcont=addcont         # 加法多项式数组
    #                                             )
    # bestfit, mpoly = uly_fit_lin_mulpol(bestfitw,                                              # 伪谱 (流量改正因子) * 权重后的 TGM 光谱
    #                                     cmp=cmp,                                               # TGM 模型字典结构
    #                                     mpoly=mpoly,                                           # 勒让德多项式字典结构
    #                                     SignalLog=SignalLog,                                   # 待测光谱字典结构
    #                                     goodPixels=goodPixels,                                 # 有效像素列表
    #                                     allow_polynomial_reduction=allow_polynomial_reduction  # 是否允许多项式阶数减少
    #                                    )

    # 6.9.2 方案 2
    # 6.9.2.1 如果 modecvg 为 None, 则设置 modecvg=0. LASP-MPFit 不使用, 但我们保留了该参数. 如有需要, 请参考 IDL 代码
    if modecvg is None:
        modecvg = 0
    # 6.9.2.2 如果 modecvg == 1 或 modecvg == 2, 则设置 itermax=499
    if (modecvg == 1) | (modecvg == 2):
        itermax = 499
    else:
        itermax = 0

    # 6.9.2.3 初始化迭代停止标志、迭代次数
    fstop, niter = 0, 0

    # 6.9.2.4 迭代循环
    while fstop == 0:
        # 6.9.2.4.1 确定每个组件的权重, LASP-MPFit 只使用了 1 个组件 (cmp)
        # 注意:
        # 1) 输入降低分辨率后的 TGM 模型光谱
        # 2) 输出低分辨率的、乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱, 更新后的 cmp、可加性多项式
        bestfitw, cmp, addcont = uly_fit_lin_weight(SignalLog,              # 待测光谱字典结构
                                                    goodPixels,             # 有效像素列表
                                                    tmp,                    # 降低分辨率后的 TGM 模型光谱
                                                    mpoly,                  # 勒让德多项式字典结构
                                                    adegree=adegree,        # 加型多项式阶数
                                                    lum_weight=lum_weight,  # 计算流量权重及其误差
                                                    cmp=cmp,                # TGM 模型字典结构
                                                    addcont=addcont         # 加法多项式数组
                                                    )
        # 6.9.2.4.2 如果 bestfitw 为 None, 则返回 bestfitw
        if bestfitw is None:
            raise ValueError('bestfitw is None!')
        # 6.9.2.4.3 weight 必须为正数, 要不然流量为非正数
        posw = np.where(cmp['weight'] > 0)[0].tolist()
        if len(posw) == 0:
            raise ValueError('All weights are not positive!')
        # 6.9.2.4.4 获取乘以伪谱 (流量改正因子) 的 TGM 光谱、更新后的勒让德多项式系数
        # 注意:
        # 1) 输入低分辨率的、乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱
        # 2) 输出低分辨率的、乘以伪谱 (流量改正因子)、权重后的 TGM 模型光谱、勒让德多项式
        bestfit, mpoly = uly_fit_lin_mulpol(bestfitw,                                              # 伪谱 (流量改正因子) * 权重后的 TGM 光谱
                                            cmp=cmp,                                               # TGM 模型字典结构
                                            mpoly=mpoly,                                           # 勒让德多项式字典结构
                                            SignalLog=SignalLog,                                   # 待测光谱字典结构
                                            goodPixels=goodPixels,                                 # 有效像素列表
                                            allow_polynomial_reduction=allow_polynomial_reduction  # 是否允许多项式阶数减少
                                            )
        # 6.9.2.4.5 归一化勒让德多项式
        # 注意: 
        # 1) 我们保留与 LASP-MPFit 一致的代码. 但下述代码块, 只是单纯用于归一化, 对参数推断结果几乎无影响 (数值精度可能对 Teff 有 0.1 K 左右影响, 其实可以删除 6.9.2.4.5)
        # 2) B(lambda_i) = c0 + c1 * P1(lambda_i) + c2 * P2(lambda_i) + ... + cn * Pn(lambda_i)
        # 3) mpoly['mpolcoefs'] 为勒让德多项式系数, 即 (c0, c1, c2, ..., cn)
        # 4) mpoly['mpolcoefs'][0] 为勒让德多项式系数的第一个系数, 即 c0
        # 5) mpoly['poly'] 为伪谱 (流量改正因子)
        # 6.9.2.4.5.1 伪谱 (流量改正因子) 除以勒让德多项式系数的第一个系数
        mpoly['poly'] /= mpoly['mpolcoefs'][0]
        # 6.9.2.4.5.2 更新 cmp 中的权重, 伪谱 (流量改正因子) 除以了勒让德多项式系数的第一个系数, 所以权重需要乘以勒让德多项式系数的第一个系数, 保证权重与伪谱 (流量改正因子) 的乘积不变
        cmp['weight'] *= mpoly['mpolcoefs'][0]
        # 6.9.2.4.5.3 勒让德多项式系数除以勒让德多项式系数的第一个系数
        mpoly['mpolcoefs'] /= mpoly['mpolcoefs'][0]

        # 6.9.2.4.5.4 LASP 使用的是 1 个 cmp, 因此跳出循环
        if len(posw) == 1:
            break

        # 6.9.2.4.5.5 更新迭代次数
        niter += 1
        # 6.9.2.4.5.6 如果、 迭代次数大于 itermax, 则停止迭代, 跳出循环
        if (abs(1.0 - (np.nansum(bestfit) / np.nansum(bestfitw))) < 5e-9) or (niter > itermax):
            fstop = 1

    # 6.10 输出伪谱 (流量改正因子)、并加权、并降低分辨率后的 TGM 光谱、勒让德多项式
    return bestfit, mpoly
