# -*- coding: utf-8 -*-
# @Time    : 2025/1/4 15:05
# @Author  : ljc
# @FileName: robust_sigma.py
# @Software: PyCharm


# 1. 简介
"""
 Python conversion of the IDL robust_sigma function .
目的:
    IDL 中定义 robust_sigma 的方法.
函数:
    1) robust_sigma
解释:
    1) robust_sigma 函数: 返回数据的鲁邦标准差.
"""


# 2. 导库
import numpy as np


# 3.鲁邦标准差
def robust_sigma(y, zero=True) -> float:

    """
        获取数据的鲁邦标准差.

        输入参数:
        -----------
        y:
          待计算标准差的数据.
        zero:
            是否考虑 0.

        输出参数:
        --------
        np.sqrt(siggma) if siggma > 0 else 0.0:
            标准差.
    """

    # 3.1 设置一个很小的数, 用于避免除以 0
    eps = 1.0e-20

    # 3.2 设置中心值
    y0 = 0.0 if zero else np.median(y)

    # 3.3 计算 MAD(median absolute deviation)
    mad = np.median(np.abs(y - y0)) / 0.6745

    # 3.4 如果 MAD 接近 0, 尝试使用平均绝对偏差
    if mad < eps:
        mad = np.mean(np.abs(y - y0)) / 0.80
    if mad < eps:
        return 0.0

    # 3.5 计算 biweight 值
    u = (y - y0) / (6.0 * mad)
    uu = u * u
    q = np.where(uu <= 1.0)[0]

    # 3.6 检查是否有足够的好点
    if len(q) < 3:
        raise ValueError('ROBUST_SIGMA: This distribution is TOO WEIRD! Returning -1!')

    # 3.7 计算总数 (处理 NaN)
    n = np.sum(np.isfinite(y))

    # 3.8 计算最终结果
    numerator = np.sum((y[q] - y0) ** 2 * (1 - uu[q]) ** 4)
    den1 = np.sum((1 - uu[q]) * (1 - 5 * uu[q]))
    siggma = n * numerator / (den1 * (den1 - 1))

    # 3.9 返回最终结果
    return np.sqrt(siggma) if siggma > 0 else 0.0


# 4. 测试
# y = np.array([1, 2, 3, 4, 0])
# result = robust_sigma(y, zero=True)  
# print(result)