chronaxis.core.moments_numpy

源码基线：8904c2f744951857659abd40aaa3535d57faf53c

1"""按有限 NumPy 块维护加权中心矩或均方根，不收集历史。23由真实有效掩码选好的值和权重按相同的二进制幂缩放。块内围绕锚点求均值和4离均差，跨块用两总体均值差合并离散度。用根均方形式存储二阶量，开方与5目标倍率在可表示范围内求值，不把两个大均方相减，也不逐样本运行标量实现。6"""78from math import frexp, hypot, ldexp, sqrt9import numpy as np10from .types import QualityFlag as Q111213class NumpyMoments:14    """持有一个锚点、相对均值、根二阶量、总权重和幂指数，跨块大小固定。"""1516    def __init__(self, operation):17        """初始化归约器。参数: 操作名称。返回: 无。异常: 调用侧校验操作。"""18        self.operation = operation19        self.weight = self.anchor = self.mean = self.root = 0.20        self.exponent = None2122    def push(self, values, weights=None):23        """整列合并活动行。参数: 有限浮点一维值，权重为正数列或等权None。返回: 无。异常: 非法权重。"""24        if not len(values):25            return26        x = np.asarray(values, dtype=np.float64)27        w = None if weights is None else np.asarray(weights, dtype=np.float64)28        if w is not None and (w.shape != x.shape or np.any(w <= 0) or not np.all(np.isfinite(w))):29            raise ValueError('二阶统计需要同长正有限权重列')30        weight = float(len(x)) if w is None else float(np.sum(w))31        peak = float(np.max(np.abs(x)))32        # 用2的整数幂共同缩放，保留高基数输入的相邻浮点差距。33        e = frexp(peak)[1] if peak else None34        if self.exponent is None or e is not None and e > self.exponent:35            if self.exponent is not None:36                shift = self.exponent - e37                self.mean, self.anchor, self.root = (ldexp(v, shift) for v in (self.mean, self.anchor, self.root))38            self.exponent = e39        exponent = self.exponent or 040        scaled = np.ldexp(x, -exponent)41        if not self.weight:42            self.anchor = float(scaled[0])43        rms = self.operation.endswith('_rms')44        if rms:45            residual, mean = scaled, 0.46        else:47            # 在全局同一锚点坐标内保存小差值；窄分布可跨块保留半ULP的均值修正。48            offsets = scaled - scaled[0]49            correction = float(np.mean(offsets) if w is None else np.sum(offsets * (w / weight)))50            mean = float(scaled[0] - self.anchor) + correction51            residual = offsets - correction52        maximum = float(np.max(np.abs(residual)))53        if maximum:54            squared = np.square(residual / maximum)55            average = float(np.mean(squared) if w is None else np.sum(squared * (w / weight)))56            block_root = maximum * sqrt(min(1., average))57        else:58            block_root = 0.59        total = self.weight + weight60        left, right = self.weight / total, weight / total61        delta = mean - self.mean62        # 平方和用hypot实现，不直接平方极小的已缩放离散度。63        combined = hypot(sqrt(left) * self.root, sqrt(right) * block_root,64                         delta * sqrt(left * right) if not rms else 0.)65        # RMS 是两块 RMS 的加权平方平均，边界夹回只修正最后一个 ULP 的舍入，66        # 保证连批常数 float64 最大值也不会越出最大输入。67        self.root = min(combined, max(self.root, block_root)) if rms else combined68        self.mean += delta * right69        self.weight = total7071    def finish(self, ddof, factor):72        """完成自由度和一次目标转换。参数: 自由度、正倍率。返回: (值或None,质量)。异常: 配置错误。"""73        if self.weight <= (ddof or 0):74            return None, Q.NONE75        if self.root == 0:76            return 0., Q.NONE77        correction = sqrt(self.weight / (self.weight - ddof)) if ddof else 1.78        mantissa, exponent = frexp(self.root * correction)79        exponent += self.exponent or 080        if self.operation.endswith('_variance'):81            mantissa, exponent = mantissa * mantissa, 2 * exponent82        f, e = frexp(factor)83        try:84            value = ldexp(mantissa * f, exponent + e)85        except OverflowError:86            return None, Q.NONFINITE87        return (value, Q.NONE) if value != 0 else (None, Q.PRECISION_LOSS)88